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1.  SUMMARY 


This  report  covers  the  work  completed  under  the  UNM/AFRL  research  grant  FA9453-1 1-1-0276 
including  the  theoretical  development  and  implementation  of  a  novel  electron  transport  method 
in  the  Geant4  toolkit,  and  other  related  work  including  the  generation  of  response  functions  for 
the  CEASE  telescope.  It  has  long  been  of  interest  to  develop  an  alternative  to  the  standard 
Condensed  History  (CH)  method  that  is  free  of  its  inherent  limitations.  The  Moment-Preserving 
method,  is  such  an  alternative  and  within  this  report  the  accuracy  and  efficiency  of  this  method 
is  demonstrated  and  contrasted  with  the  CH  method.  The  Moment-Preserving  (MP)  Monte  Carlo 
electron  transport  method  was  developed  and  implemented  in  the  Geant4  toolkit  as  C++  classes 
for  constructing  the  libraries  that  are  required  by  the  MP  method  physics  models.  Given  the 
successful  implementation  of  the  physics  models,  the  MP  method  was  then  tested  on  a  wide 
variety  of  problems  including  a  theoretical  test  suite,  experimental  validations,  and  the  CEASE 
particle  telescope. 


2.  INTRODUCTION 

This  report  covers  the  work  completed  under  the  UNM/AFRL  research  grant  FA9453-1 1-1-0276 
during  the  period  of  April  2011  to  March  2015  including  the  primary  objective,  or  the  theoretical 
development  and  implementation  of  a  novel  electron  transport  method  in  the  Geant4  toolkit,  and 
other  related  work  including  the  modeling  of  detectors  associated  with  the  AFRL  Demonstration 
and  Science  Experiment  (DsX)  detectors  for  use  in  response  function  determination.  In  efforts  to 
complete  the  primary  objective,  the  Moment-Preserving  (MP)  Monte  Carlo  electron  transport 
method  was  developed  and  implemented  in  the  Geant4  toolkit.  Specifically,  C++  classes  were 
developed  for  constructing  the  reduced  order  physics  (ROP)  differential  cross-section  (DCS) 
libraries  that  are  required  by  the  MP  method  physics  models  developed  for  the  Geant4  toolkit. 
Given  the  successful  implementation  of  the  physics  models,  the  MP  method  was  then  tested  on  a 
wide  variety  of  problems  including  a  theoretical  test  suite,  experimental  validations,  and 
generation  of  response  functions  for  the  CEASE  telescope  (i.e.  a  DsX  instrument).  The 
motivation  for  MP  method  and  the  key  properties  are  now  summarized. 

Analog  electron  transport  is  computationally  impractical;  and  therefore,  methods  such  as 
condensed  history  were  developed  to  alleviate  the  computational  cost  associated  with  analog 
Monte  Carlo  electron  transport.  However,  the  condensed  history  method  suffers  from  limitations 
inherent  to  the  method  resulting  from  the  formulation  of  the  method.  Therefore,  it  is  of  interest  to 
develop  an  advanced  method  that  is  free  of  the  limitations  of  condensed  history,  while  remaining 
efficient  and  accurate  with  respect  to  the  analog  transport  model.  The  analog  description  of 
transport  can  be  mathematically  expressed  in  terms  of  the  linear  Boltzmann  or  transport  equation 
for  the  phase  space  particle  distribution  function,  where  interaction  physics  is  represented 
through  total  cross  sections  and  differential  cross  sections  for  energy  transfer  and  angular 
deflection.  In  the  so-called  advanced  method,  a  reduced  order  transport  problem  is  posed,  but 
with  modified  cross  sections  such  that  the  resulting  single-event  Monte  Carlo  simulation  is 
computationally  efficient  (minutes  vs.  days).  Underlying  the  approach  is  a  method  of 
systematically  preserving  energy-loss  and  angular  deflection  moments  of  the  analog  differential 


Approved  for  public  release;  distribution  is  unlimited. 

1 


cross  sections  and  using  only  a  finite  number  of  these  moments  to  reconstruct  a  reduced  order 
physics  differential  cross  section.  This  method  is  characterized  by  the  following  properties: 

1 .  The  mean  free  path  of  the  reduced  order  physics  model  will  be  longer  than  the  true  or 
analog  mean  free  path  and  the  differential  cross  sections  less  peaked,  with  accuracy 
depending  on  the  number  of  moments  preserved. 

2.  The  moment-preserving  methodology  will  enable  accuracy  to  systematically  approach 
analog  by  incorporating  increasingly  higher  order  moments  with  precision. 

3.  By  retaining  the  “look  and  feel”  of  the  analog  transport  process,  i.e.,  using 
approximations  based  on  Boltzmann  or  integral  fonns  of  the  collision  operators 
describing  interactions,  collisions  are  exactly  treated  as  Markov  events,  thereby  allowing 
material  and  vacuum  interfaces  to  be  handled  without  algorithmic  modifications  in 
condensed  history  methods. 

4.  In  aggregate,  the  smoother  cross  section  variation  with  energy  and  angle  as  well  as  the 
correct  transport  mechanics  will  result  in  considerable  computational  savings  over  using 
strictly  analog  Monte  Carlo  simulation.  We  expect  the  savings  to  be  such  that  the 
moment-preserving  method  will  be  computationally  competitive  with  condensed  history, 
but  it  also  has  the  potential  to  be  more  efficient  than  condensed  history. 

5.  The  methodology  will  be  independent  of  charged  particle  species  and  will  therefore 
apply  to  species  of  disparate  masses,  such  as  protons  and  electrons. 

6.  The  methodologies  proposed  here  can  be  directly  incorporated  into  production  Monte 
Carlo  codes  such  as  MCNPX  and  Geant4  without  necessitating  new  logic  because  this 
approach  will  treat  charged  particles  similar  to  neutral  particles  with  regards  to  transport. 

Given  a  summary  of  the  advanced  Monte  Carlo  method  and  its’  properties,  the  background 
information  necessary  to  understand  the  motivation  of  developing  such  a  method  is  provided. 

3.  BACKGROUND 

The  need  for  computational  charged  particle  transport  developed  from  early  work  in  charged 
particle  transport  theory,  which  emerged  as  a  flourishing  sub-branch  of  mathematical  physics 
when  fast  charged  particles  became  available  to  the  experimentalist  [2].  As  computer 
technology  improved,  the  problems  of  interest  to  charged  particle  computational  physicist 
expanded  to  areas  including:  space  physics,  accelerator  physics,  medical  physics,  health  physics, 
and  electromagnetic  pulses.  The  advantage  of  computational  charged  particle  transport  over 
analytical  transport  is  the  possibility  of  simulating  complicated  geometries  and  sophisticated 
boundary  conditions  or  source  configurations,  which  are  all  characteristics  of  real  world 
applications.  In  other  words,  it  is  possible  to  simulate  real,  physical  phenomena  using  charged 
particle  transport  codes.  An  example  of  such  a  code  is  the  Geant4  toolkit  [3]  which  is  used 
frequently  on  problems  including:  design  of  full-scale  experiments  such  as  the  Large  Hadron 
Collider  [4,  5],  design  of  radiation  therapy  machines  [6]  as  well  as  treatment  planning  systems 
[7],  estimation  of  detector  geometric  factors  [8],  shielding  calculations  [9],  and  EMP 
calculations  [10],  It  is  undisputed  that  particle  transport  codes  play  an  important  role  in  the 
research  and  development  of  charged  particle  applications  and  reasonable  to  suggest  that 
particle  transport  codes  will  continue  to  play  an  important  role  into  future.  Therefore, 
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algorithmic  improvements  to  particle  transport  codes  are  critical  to  improving  the  field  of 
computational  charged  particle  transport. 

To  make  clear  the  impact  of  this  work  and  the  importance  of  the  remaining  chapters,  it  is 
necessary  to  provide,  at  least,  a  cursory  discussion  on  what  is  meant  by  charged  particle 
transport  codes  and  the  associated  challenges.  First,  the  purpose  of  charged  particle  transport 
codes  is  to  obtain  solutions  to  the  Boltzmann  trans-  port  equation  [11]  using  stochastic 
methods  like  Monte  Carlo  [12]  or  deterministic  methods  like  Sjy  [13].  The  Boltzmann  transport 
equation  is  a  balance  equation  for  particles  in  a  six-dimensional  phase-space  including  space, 
angle,  and  energy.  The  solution  to  this  equation  describes  the  particle  population  and  is  referred 
to  as  the  angular  flux.  For  a  given  analog  DCS  model,  the  corresponding  transport  equation  is 
referred  to  as  the  analog  model  and  the  angular  flux  is  assumed  to  be  exact.  Analog,  (detailed, 
step  by  step),  simulation  is  feasible  under  strict  circumstances  (relatively  low  energies,  thin 
targets,  ...),  but  for  high-energy  electrons  (above  a  few  hundred  keV),  the  number  of 
interactions  suffered  by  an  electron  along  its  trajectory  is  too  large  for  detailed  simulations 
[14].  For  this  reason,  numerous  computationally  efficient  approximate  methods  have  emerged 
over  the  past  60  years.  The  most  notable  and  prolific  approximate  method  is  Condensed  History 
(CH).  Berger  describes  CH  as  an  artificially  constructed  random  walk,  each  step  of  which  takes 
into  account  the  combined  effects  of  many  collisions  [2].  The  distances  between  collisions  or  the 
steps  are  significantly  longer  than  those  associated  with  the  analog  problem,  making  CH 
efficient.  However,  the  theoretical  basis  and  practical  implementation  of  the  CH  algorithm 
introduces  inherent  and  irreducible  limitations  that  are  unique  to  the  method. 

It  is  of  interest  to  develop  a  method  free  of  such  limitations.  To  better  understand  benefits  of 
such  a  method,  the  remaining  sections  provide  a  qualitative  discussion  of  the  literature  relevant 
to  the  analog  problem  and  the  associated  analog  differential  cross-sections  (DCS),  the 
Condensed  History  (CH)  method,  and  Reduced  Order  Physics  (ROP)  models  including  the 
Moment-Preserving  (MP)  method.  The  analog  problem  is  the  point  of  departure  for  both  CH 
and  ROP  models,  so  we  begin  with  a  discussion  of  the  analog  problem. 

3.1.  The  Analog  Problem 

The  analog  description  of  transport  can  be  mathematically  expressed  in  terms  of  the  linear 
Boltzmann  or  transport  equation  for  the  angular  flux,  where  the  interaction  physics  are 
represented  through  total  cross  sections  (inverse  mfp)  and  the  DCSs  for  angular  deflection 
and  energy  loss.  The  total  cross  sections  and  DCSs  appear  in  the  elastic  and  inelastic  collision 
operators,  which  are  integral  operators  or  Boltzmann-type  operators.  Though  electrons  can 
undergo  several  different  electromagnetic  interactions,  the  dominant  interactions  are  elastic 
collisions  with  nuclei  and  inelastic  collisions  with  atomic  electrons.  Typical  DCSs  for  elastic 
nuclear  scattering  include  relativistic  screened  Rutherford,  Wentzel,  or  the  partial-wave 
expansion  [15,  16,  17],  while  typical  DCSs  for  inelastic  electronic  scattering  include 
Rutherford,  Moller,  or  the  Evaluated  Electron  Data  Library  [18,  19].  These  DCSs  are  highly 
peaked  about  small  changes  in  direction  and  small  energy  losses  and  the  associated  total 
cross  sections  are  very  large  resulting  in  extremely  short  mfps.  Interaction  physics  of  this 
nature  present  a  difficult  computational  task. 
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Boundary  conditions  for  the  transport  equation  depend  on  the  application,  but  typically  include 
vacuum,  pencil  beams,  or  sources  distributed  in  space,  energy,  and  angle.  However,  mono- 
energetic  pencil  beams  are  very  common  in  electron  trans-  port  and  further  complicate  the 
computational  challenges  because  pencil  beams  are  singular  in  space,  energy,  and  angle. 

The  problem  of  computational  inefficiencies  associated  with  the  analog  physics  was 
recognized  immediately  by  early  charged  particle  computational  physicists.  In  fact,  Berger  [2] 
acknowledged  that  a  direct  analog  Monte  Carlo  procedure  would  be  quite  costly,  because  of  the 
enormous  number  of  collisions  that  must  be  sampled.  For  example,  it  takes  on  the  order  of  tens 
of  thousands  of  collisions  for  an  electron  with  energy  of  roughly  1-MeV  to  slow  down  to  1-keV, 
while  only  20  to  30  Compton  scatterings  will  reduce  the  energy  of  a  photon  from  several  MeV 
down  to  1-keV  or  18  elastic  collisions  in  hydrogen  will  reduce  a  neutron  from  2-MeV  to 
thermal  energies  [2],  In  1963,  only  one  calculation  by  direct  analog  Monte  Carlo  was  reported 
[20],  Since  then,  several  analog  Monte  Carlo  electron  transport  codes  have  been  developed 
[14,  21]  and  a  few  production  codes  have  included  analog  physics  options  [22,  23,  24]. 

Solutions  to  the  analog  problem  are  exact  for  a  given  DCS;  therefore,  it  is  understandable  that 
analog  physics  options  were  implemented,  despite  the  fact  that  analog  Monte  Carlo  is 
computationally  inefficient.  Moreover,  it  is  feasible  to  use  analog  Monte  Carlo  for  occasional 
calculations  if  significant  computing  resources  are  available.  However,  analog  Monte  Carlo 
remains  impractical  for  routine  calculations  with  exception  of  very  restrictive  problems  like 
transport  through  optically  thin  materials.  For  this  reason,  approximate  methods  remain  a 
critical  component  of  most  Monte  Carlo  electron  transport  codes. 

3.2.  Condensed  History 

Condensed  history  has  been  the  prevailing  approximate  method  in  computational  charged 
particle  transport  since  the  emergence  of  the  field.  The  usual  practice  is  to  use  "condensed" 
(class  I)  simulation  methods,  in  which  the  global  effect  of  multiple  interactions  is  described  by 
means  of  approximate  multiple  scattering  theories.  Alternatively,  one  can  use  "mixed"  (class  II) 
schemes  in  which  hard  (catastrophic)  interactions,  with  energy  loss  or  angular  deflection  above 
given  thresholds,  are  simulated  individually.  For  a  given  set  of  DCSs,  class  II  schemes  are 
intrinsically  more  accurate  than  class  I  simulations  [14].  Some  examples  of  codes  containing 
class  I  schemes  are  ETRAN,  ITS,  MCNP,  while  examples  of  production  codes  containing  class 
II  CH  schemes  are  Geant4,  PENELOPE,  and  EGS4  [3,  24,  25],  Both  class  I  and  class  II 
schemes  utilize  various  results  from  multiple  scattering  theory,  which  is  a  subbranch  of 
mathematical  physics  developed  around  the  solution  of  the  trans-  port  equation  with  limited 
applicability  resulting  from  severe  restrictions  required  to  obtain  analytical  solutions  [26,  27, 

28,  29].  The  analytical  solutions  or  multiple  scattering  distributions  describe  the  angular  or 
energy  distributions  of  electrons  after  traveling  some  distance  s  or  a  step,  that  are  on  the  order  of 
hundreds  of  analog  mfps. 

The  major  distinction  between  class  I  and  class  II  schemes  is  how  the  grouping  of  collisions  is 
handled.  That  is,  class  I  schemes  utilize  precomputed  multiple-scattering  (MS)  distributions 
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[30,  31,  32,  33]  determined  for  fixed  step-sizes  on  a  fixed  energy  grid.  For  this  reason,  energy 
straggling  is  sampled  from  a  MS  distribution  [32,  33]  and  secondaries  are  accounted  for  on  an 
average  basis.  Therefore,  it  is  not  possible  to  distinguish  between  inelastic  collisions  resulting  in 
secondary  production  (hard)  and  those  that  do  not  (soft).  In  contrast,  a  class  II  scheme  like  EGS 
allows  all  physical  processes  and  boundaries  to  affect  the  choice  of  step  size  [34].  Thus,  distance 
to  hard  collision  is  exponentially  distributed,  while  secondary  production  is  treated  exactly  by 
sampling  energy-loss  from  the  inelastic  DCS  above  the  secondary  production  threshold. 

Regardless  of  the  choice  in  the  scheme,  substantial  efficiency  gains  over  analog  Monte  Carlo 
can  be  realized  with  CH.  However,  the  accuracy  of  early  fonns  of  CH  was  strongly  dependent 
on  step-size  and  while  it  was  found  that  reducing  the  electron  step-size  causes  the  results  to 
converge  to  the  correct  values,  the  computing  time  increases  rapidly  in  proportion  to  the  inverse 
of  the  step-size  [35,  36,  37].  Therefore,  special  algorithms  like  PRESTA  [38]  were  developed  to 
select  the  optimal  step-size  during  the  process  of  a  Monte  Carlo  simulation.  Without  an 
algorithm  like  PRESTA,  one  must  resort  to  a  tedious  study  to  detennine  the  optimal  step-size 
such  that  acceptable  accuracy  and  efficiency  is  achieved.  However,  this  optimization  may  not  be 
universal.  The  various  production  codes  currently  available  differ  in  this  optimization  issue. 
Some  codes,  like  ITS,  have  pre-detennined  step-size  parameters,  while  codes  like  EGS  utilize 
the  PRESTA  algorithm  [38]  or  random  hinging  combined  with  lateral  corrections  found  in 
PENELOPE.  In  addition  to  step-size  limitations,  condensed  history  suffers  from  inconsistent 
handling  of  the  material  and  free  surface  boundaries.  Material  interfaces  are  a  fundamental 
challenge  for  condensed  history  because  the  MS  distributions  are  infinite  medium  solutions  and 
are  not  valid  for  heterogeneous  regions.  Therefore,  if  a  material  interface  is  encountered,  a 
special  algorithm  like  the  Jordan-Mack  algorithm  [39]  or  PRESTA  is  required.  Another  issue 
specific  to  class  I  schemes  that  utilize  the  Goudsmit-Saunderson  distribution  [30,  31]  for 
sampling  angular  deflection  is  that  the  numerical  methods  required  to  generate  the  Goudsmit- 
Saunderson  distribution  are  sensitive  to  small  step-sizes.  The  backwards  recurrence  used  to 
generate  the  Goudsmit-Saunderson  is  unstable  for  small  steps.  Even  if  the  step-size  is 
sufficiently  long  that  the  backwards  recurrence  is  stable,  it  is  possible  that  more  than  the  pre¬ 
scripted  number  of  recurrence  coefficients  are  required  to  accurately  resolve  the  Goudsmit- 
Saunderson  distribution  [40] . 

3.3.  Reduced  Order  Physics  Models 

For  the  reasons  stated  above,  various  alternatives  to  the  CH  method  were  developed.  Of 
particular  interest  are  alternatives  referred  to  as  Reduced  Order  Physics  (ROP)  models,  which 
are  a  family  of  transport-based  approximations.  ROP  models  include  various  approximate 
representations  of  the  analog  collision  operators  that  are  of  both  integral  and  differential  fonns. 
ROP  models  are  obtained  through  some  type  of  regularization  procedure  that  removes  reduce 
the  nearly-singular  behavior  of  the  analog  DCS,  while  systematically  capturing  the  key  physics 
through  preservation  of  angular  and  energy-loss  moments  of  the  analog  DCSs.  Moreover,  the 
zeroth  or  the  total  cross  section  is  not  preserved,  rather  it  is  detennined  self-consistently  by  the 
method.  The  resulting  ROP  model  is  then  characterized  by  less  peaked  scattering  with  mfps 
longer  than  the  analog  mfp.  Efficiency  is  achieved  by  not  preserving  the  total  cross  section, 
while  accuracy  is  achieved  by  preserving  the  necessary  number  of  moments  beyond  the  zeroth. 
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There  are  numerous  approximations  that  qualify  as  ROP  models  including  Fokker-Planck, 
Boltzmann  Fokker-Planck,  Generalized  Fokker-Planck,  and  Generalized  Boltzmann  Fokker- 
Planck.  The  remainder  of  this  chapter  is  devoted  to  introducing  each  of  these  approximations 
and  indicating  how  these  methods  contributed  to  the  development  of  the  MP  method. 

The  classical  Fokker-Planck  operator  is  obtained  by  Taylor  expanding  the  scattering  kernels. 
This  is  a  reasonable  approach,  assuming  the  angular  flux  is  sufficiently  smooth,  because  the 
DCSs  fall  of  rapidly  away  from  small  deflection  cosines  and  energy-losses.  The  resulting 
operator  is  differential  in  angle  and  energy  and  models  elastic  and  inelastic  scattering  as 
diffusive  processes.  As  a  result,  the  Fokker-Planck  operator  does  not  capture  large-angle  scatter 
and  can  lead  to  energy  gains.  Pomraning  confirmed  these  inconsistencies  by  showing  that  the 
Fokker-Planck  operator  is  an  asymptotic  limit  of  the  Boltzmann  collision  operator  and  only 
valid  for  unrealistically  peaked  scattering  [41].  That  is,  the  FP  approximation  is  strictly  valid  in 
the  limit  that  the  total  scattering  cross  section  goes  to  infinity  and  the  mean  deflection  cosine 
goes  to  unity  such  that  the  transport  cross  section  remains  bounded  [42].  Under  these  conditions, 
large  angle  scattering  is  negligible  and  the  FP  operator  is  equivalent  to  the  Boltzmann  integral 
collision  operator.  Clearly,  there  is  a  limitation  on  the  type  of  physics  for  which  the  classical  FP 
operator  is  valid.  Regardless,  various  implementations  of  Fokker-Planck  operator  were  studied 
for  in  detenninistic  settings  [43,  44,  45,  46]. 

In  efforts  to  incorporate  large-angle  scattering,  a  kernel  decomposition  approach  was  introduced 
by  Ligou  and  is  referred  to  as  the  Boltzmann-Fokker-Planck  (BFP)  equation  [47,  48,  49].  Ligou 
recognized  that  it  is  easier  to  numerically  treat  forward-  peaked  elastic  scattering  and  small 
energy  losses  associated  with  such  scattering  using  Fokker-Planck  (FP)  differential  operators 
[45]  rather  than  Boltzmann  integral  operators.  However,  the  FP  operator  does  not  accurately 
capture  large  angle  scattering.  Therefore,  it  is  necessary  to  decompose  the  scattering  cross 
section  into  its  singular  and  smooth  components  and  apply  the  FP  approximation  to  the 
singular  component  while  leaving  the  smooth  operator  intact  [50].  One  important  feature  of  the 
decomposition  process  is  that  there  is  no  rigorous  definition  of  the  components  so  there  are 
infinite  decompositions.  The  key,  as  indicated  by  Landesman  and  Morel  [50],  is  to  select  a 
decomposition  that  is  not  only  accurate  and  efficient,  but  also  easily  integrated  into  existing 
transport  codes.  The  early  methods  for  solving  the  BFP  equation  were  deterministic,  but  Morel 
and  Sloan  [51,  52]  developed  a  hybrid  multigroup/continuous-energy  Monte  Carlo  method  for 
solving  the  BFP  equation  or  the  MGBFP  method.  Morel  described  the  MGBFP  method  as  a 
new  form  of  condensed  history  with  the  major  distinction  being  that  path  lengths  between 
collision  sites  are  exponentially  distributed. 

In  retrospect,  application  of  kernel  decomposition  serves  to  stabilize  the  divergent  behavior  of 
the  Fokker-Planck  expansion  [53]  by  effectively  renonnalizing  the  expansion  coefficients.  This 
renormalization  process  is  central  to  the  Generalized  Fokker-Planck  (GFP)  method  [54,  55]  and 
resulted  after  attempting  to  generalize  the  FP  operator  to  a  broad  class  of  physics  with  higher 
order  FP  operators  [53].  In  principle,  the  FP  approximation  can  be  improved  by  retaining 
higher  order  tenns  in  the  Taylor  expansion  of  the  collision  operators.  Thus,  leading  to  a  more 
accurate  description  of  large  angle  scattering,  but  this  is  only  true  for  specific  kernels.  For 
example,  it  is  known  that  there  are  no  valid  FP  operators  for  the  Henyey-Greenstein  kernel  [41] 


Approved  for  public  release;  distribution  is  unlimited. 
6 


and  that  the  standard  FP  operator  is  only  marginally  valid  for  screened  Rutherford  or  screened 
Mott  [42],  However,  Pomraning  [53]  completed  the  same  kernel  analysis  for  the  exponential 
and  delta  function  and  showed  that  higher  order  FP  operators  are  valid  for  kernels  such  as  these, 
but  they  are  unphysical.  In  addition  to  characterizing  the  stability  of  higher  order  FP  operators 
with  respect  to  specific  scattering  kernels,  this  analysis  revealed  the  explicit  role  of  the  analog 
DCS  moments  to  developing  approximations.  That  is,  higher  order  FP  expansions  suggest  that 
the  accuracy  of  the  approximation  is  improved  by  incorporating  higher  order  moments.  GFP 
approaches  provide  a  means  to  stabilize  higher-order  FP  operators  through  renormalization  of 
the  various  terms  that  appear  in  the  operator.  The  renormalization  process,  in  essence,  allows 
for  an  arbitrary  ordered  FP  operator;  hence,  the  approach  referred  to  as  generalized  FP.  Lastly, 
work  on  GFP  showed  that  the  key  to  a  stable  ROP  model  is  preserving  the  integral  fonn  of  the 
collision  operator  and  retaining  an  arbitrary  number  of  low  order  moments  while  approximating 
all  of  the  higher  order  moments. 

Insight  gained  from  the  foregoing  approximations  -  FP,  BFP,  and  GFP  -  combined  with  Lewis 
theory  [56,  57]  and  practical  experience  from  implementing  and  testing  CH  led  to  the 
development  of  the  Generalized  Boltzmann  Fokker-Planck  (GBFP)  method  [58],  Like  the  other 
ROP  models,  the  GBFP  method  is  a  transport-based  approximation.  However,  the  GBFP 
method  is  more  appropriate  for  Monte  Carlo  calculations  like  CH.  The  emphasis  of  GBFP 
method  is  the  development  of  stable,  moment-preserving  representations  of  the  collision 
operators.  This  is  achieved  by  constructing  ROP  DCSs  that  are  moment-preserving,  per  Lewis 
theory,  while  leaving  the  collision  operators  intact  by  simply  replacing  the  analog  DCSs  with 
the  ROP  DCSs.  Another  key  feature  of  the  GBFP  approach  is  that  angular  deflection  and 
energy-loss  interactions  depend  only  on  Legendre  moments  and  energy-loss  moments,  not  the 
detailed  form  of  the  analog  DCSs.  Therefore,  the  GBFP  method  is  applicable  to  both  continuous 
DCSs  [59]  and  tabulated  DCS  data  [60]  and  code  modifications  are  not  required  for  physics 
refinements  because  this  is  accomplished  by  providing  moments  of  the  desirable  analog  model 
through  input  data.  The  bulk  of  the  work  on  the  GBFP  method  to  date  emphasizes  the  most 
simple,  but  also  the  most  effective  ROP  DCSs  referred  to  as  the  discrete  and  hybrid  models. 

The  discrete  model  [61]  is  a  superposition  of  discrete  points  and  weights  over  the  full  range  of 
the  DCS,  while  the  hybrid  model  [62]  is  a  superposition  of  discrete  points  and  weights  over  the 
peaked  portion  of  the  DCS  and  the  tail  is  represented  exactly  by  the  analog  DCS.  Recently,  the 
name  Moment-Preserving  method  was  adopted  in  place  of  the  GBFP  method  to  distinguish  the 
method  from  Fokker-Planck  and  to  emphasize  the  concept  most  fundamental  to  the  method  - 
moment  preservation. 

4.  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

In  this  section,  the  Moment-Preserving  method  is  presented  along  with  a  discussion  of  the  analog 
DCSs  and  the  analog  transport  model. 
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4.1.  Analog  Differential  Cross-Sections  and  the  Analog  Transport  Model 

In  the  following  sections,  the  analog  physics  are  presented  along  with  the  analog  transport 
model.  The  analog  electron  physics  render  application  of  the  Monte  Carlo  method  to  analog 
transport  model  computationally  inefficient.  Specific  features  associated  with  the  analog  physics 
that  attribute  to  the  extreme  computational  cost  are  discussed.  Herein,  solutions  to  the  analog 
transport  model  are  referred  to  as  the  analog  benchmark  and  accuracy  and  efficiency  of  the 
Moment-Preserving  method  are  measured  with  respect  to  the  analog  benchmark.  Therefore,  the 
analog  transport  model  is  presented  along  with  boundary  conditions  and  applicable  assumptions. 

4.1.1.  Elastic  and  Inelastic  Analog  Differential  Cross-Sections 

In  this  paper,  we  are  concerned  with  the  interactions  that  render  analog  Monte  Carlo  electron 
transport  computationally  inefficient.  That  is,  elastic  collisions  with  atomic  nuclei  (primary 
source  of  deflection)  and  inelastic  collisions  with  atomic  electrons  (primary  source  of  energy- 
loss).  These  interactions  are  characterized  by  highly  peaked  DCSs  about  small  changes  in 
direction  and  small  energy  transfers  or  energy  losses  resulting  in  extremely  large  total  cross- 
sections  or  short  mfps  on  the  order  of  microns. 


Figure  1:  Electron  interaction  diagrams  for  elastic  and  inelastic  scattering  [24]. 

The  analog  DCSs  used  in  this  work  are  referred  to  as  the  partial-wave  elastic  DCS  and  the 
Moller  inelastic  DCS.  However,  the  nature  of  electronic  interactions  is  understood  through 
simple  forms  of  the  elastic  and  inelastic  DCSs  such  as  screened  Rutherford  and  the  Rutherford 
energy-loss  DCS.  These  DCSs  describe  Coulombic  interactions  between  incident  and  target 
particles  and  are  therefore  dominated  by  distant  collisions.  As  such,  the  DCSs  are  highly  peaked 
and  near-singular  about  forward  directions  and  zero  energy-losses  resulting  in  extremely  short 
collision  mean  free  paths.  For  instance,  nuclear  scattering  of  high  energy  electrons  is  well 
described  by  the  screened  Rutherford  DCS,  or 


=  2nr2  Z2N 


1  -/?2  1 

/?4  [1  +  2  ry-Mo]2' 


(1) 
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where,  aside  from  physical  parameters,  /?  —  v/c,v  being  the  speed  of  the  ion  and  c  the  speed  of 
light  in  vacuum,  and  r]  is  the  screening  parameter  typically  in  the  range  of  10-7  to  10~4  with 
smaller  numbers  corresponding  to  higher  energies  and  lighter  target  nuclei.  From  this  it  follows 
that  the  total  scattering  cross  section 


ZelW 


17(F) 


»  1, 


demonstrating  that  the  mean  free  path  is  very  small,  while 


Eei(E,p  o  =  !)  1 

Iei(F,Mo<0.95)~  772(£)  ' 


(2) 

(3) 


Demonstrating  that  scattering  is  highly  forward  peaked.  Thus,  the  picture  of  charged  particle 
interactions  is  one  of  infinitely  frequent  collisions  with  infinitesimally  small  changes  of  state. 
Clearly  then,  except  in  rare  cases  of  large  angle  scattering,  a  single  collision  is  not  representative 
of  the  transport  process.  Noticeable,  indeed  physically  interesting,  changes  in  the  particle  state 
occur  on  the  average  only  after  a  large  number  of  collisions.  A  more  useful  measure  of  the  rate 
of  change  of  the  particle  state  is  then  provided  by  the  momentum  transfer  or  transport  cross 
section  defined  by 


StrOO 

from  which  it  readily  follows  that 


=  2x[  a- 


Mo)^ei(E,  p0)dp0, 


gtrCg) 

Zel(E) 


~r](E)ln 


\  1  1 

im. 


«  1. 


(4) 

(5) 


The  physical  significance  of  this  result  is  that  it  is  on  the  scale  of  a  transport  mean  free  path  that 
an  initially  collimated  angular  distribution  will  have  broadened  considerably.  Key  to  any 
approximate  treatment  of  charged  particle  transport  is  then  to  recognize  that  it  is  not  necessary  to 
resolve  the  solution  with  respect  to  the  true  mean  free  path.  More  realistic  cross  sections  such  as 
the  PW  DCSs  differ  from  the  screened  Rutherford  cross  section  at  larger  scattering  angles,  but 
the  singular  behavior  at  forward  scattering  angles  is  given  precisely  by  screened  Rutherford. 
Likewise,  proton  and  heavier  ion  scattering  is  described  by  the  screened  Rutherford  DCS  near 
forward  scattering  angles  and  all  of  the  above  considerations  apply. 


Inelastic  energy  losses  are  similarly  singular  about  zero  energy  transfers.  For  instance,  the 
energy-loss  DCS  for  an  ion  colliding  with  a  target  electron  is  accurately  described  by  the 
relativistic  Rutherford  cross  section  given  by 


Zln(E,Q) 


0.1536 


Z\Z2p  1 
a2  p2Q2 


Q_ 

I 

max 


min 


<  Q  <  Qmax(P), 


(6) 


Approved  for  public  release;  distribution  is  unlimited. 

9 


where  Q  is  the  single  collision  energy  loss  variable.  We  see  from  this  that  the  total  cross  section 
is 


Zfn(E) 


-f 


Zin(E,Q)dQ 


P2Qr 


while  mean  energy  loss  per  unit  pathlength  travelled,  or  the  stopping  power,  is  given  by 


S(E) 


Jr 


Qmax 

Qmin 


Q  Zin(E,Q)dQ  ~  In 


m 


(7) 

(8) 


Thus,  because  Qmtn  «  Qmax>  we  find  for  fixed  /?  that  Ein(E)  »  S(E)  and  hence  the  mean 

5(£) 

energy  loss  per  collision  (AE)  =  — —  «  1.  As  an  example,  for  protons  incident  on  tungsten, 

Sin(E) 

the  total  interaction  cross  section  ranges  from  about  2  x  103  cm~ 1  at  an  energy  of  2  GeV  to 
about  8  x  104  cm~ 1  at  10  MeV,  while  the  mean  energy  loss  varies  over  the  same  energy  range 
from  6  x  10~3  MeV  to  2.5  x  10-3  MeV.  Given  that  the  range  of  the  proton  is  about  30  cm,  the 
difficulty  of  resolving  energy  spectra  or  dose  distributions  using  an  analog  description  of  the 
physics  can  be  readily  appreciated.  Electron  and  positron  inelastic  electronic  collisions  are 
similarly  described  by  DCSs  that  are  singular  about  zero  energy  transfers,  specifically  given  by 
Moller  and  Bhabha  cross  sections  respectively  for  electrons  and  positrons. 


That  is  to  say,  the  physical  interactions  that  dominate  scattering  and  energy-loss  of  all  charged 
particles  are  mediated  by  long  range  Coulomb  forces  which  strongly  weight  small  deflections 
and  small  energy-losses  per  collision.  As  a  consequence,  solving  the  analog  problem  by  Monte 
Carlo  or  deterministic  methods  is  completely  impractical.  Significant  changes  to  the  particle 
distribution  function  in  space,  angle,  and  energy  can  arise  only  from  the  cumulative  effect  of  a 
large  number  of  collisions  and  any  approximate  transport  description  must  recognize  and  exploit 
this  fact. 


In  this  work,  elastic  and  inelastic  scattering  are  given  by  the  partial-wave  DCS  and  the  Moller 
DCS  respectively.  The  PW  DCSs  are  numerically  evaluated  using  the  ELSEPA  code  [17]  and 
are  considered  accurate  representations  of  the  elastic  scattering  of  electrons  by  atomic  nuclei  at 
energies  above  roughly  1-keV.  The  PW  DCSs  utilized  in  this  work  are  similar  to  those  included 
in  the  ICRU-ELSEP  database  [63]  with  an  exception  being  that  the  PW  DCSs  were  evaluated  on 
an  energy  grid  with  16  logarithmically  spaced  points  between  10-MeV  and  20-MeV  and  107 
logarithmically  spaced  points  between  1-keV  and  10-MeV.  A  log-log  linear  interpolation  scheme 
[64]  was  used  to  obtain  both  the  total  cross  section  and  the  DCS  for  a  given  energy  between  grid 
points. 

The  elastic  collisions  with  nuclei  considered  herein  deflect  the  incident  electron  through  some 
scattering  angle.  While  some  energy  is  transferred  to  the  nucleus  during  an  elastic  collision,  the 
mass  of  the  nucleus  is  so  large  in  comparison  to  the  electron  mass  that  the  energy  transferred  to 
the  nucleus  is  assumed  negligible.  As  mentioned,  elastic  scattering  DCSs  are  extremely  peaked 
about  very  small  deflection  angles  and  in  some  cases  vary  up  to  28  orders  of  magnitude  resulting 
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in  extremely  short  mfps  as  small  as  fractions  of  a  micron.  This  behavior  is  seen  in  Figure  2 
where  PW  DCSs  are  presented  for  aluminum  and  gold  at  various  energies. 


0.0001  0.01  1  10  50  100  150  0.0001  0.01  1  10  50  100  150 


Figure  2:  Partial-wave  elastic  differential  cross-sections  for  various  energy  electrons 
scattering  with  aluminum  (left)  and  gold  (right)  nuclei. 


The  Moller  inelastic  DCS  [18]  accounts  for  energy  transferred  to  atomic  electrons  resulting  in 
ionization  or  excitation  of  the  target  atom.  In  the  Moller  DCS,  binding  energies  for  the  electrons 
shells  are  neglected;  therefore,  Moller  is  an  approximation  over  all  of  the  shells  where  the 
minimum  energy  transfer  Qmtn  is  assumed  to  be  the  mean  ionization  potential  iMev  [65].  The 
resulting  inelastic  DCS  in  which  the  incident  electron  with  rest  mass  of  and  kinetic  energy  E 
transfers  energy  Q  to  the  slower  electron  is  given  by 


ZinV'Q) 


2nrezZm0c2 


'  1  1  _  1 
Q2  +  {E  -  QY  +  (E  +  m0c2) 


m0c2(2E  +  m0c2) 
Q(E-QXE  +  m0c2)  \’ 


(9) 


where  Q  6  [ Qmin ,  Qmax]  and  the  upper  bound  is  Qmax  =  E/2.  Like  the  elastic  scattering  DCS, 
the  inelastic  DCS  is  peaked  about  small  energy  transfers  or  Qmtn-  The  Moller  DCS  is  not 
considered  an  accurate  representation  of  inelastic  collisions  except  for  large  energy  transfers.  In 
fact,  the  first  moment  of  the  Moller  DCS,  that  is  the  stopping  power,  fails  to  reproduce  the 
stopping  powers  given  by  the  ICR  Report  37  [66].  Here,  we  renonnalize  the  Moller  DCS  such 
that  the  first  moment  is  in  agreement  with  the  stopping  powers  in  ICRU  Report  37. 


4.1.2.  Analog  Transport  Model 

Given  a  description  of  the  electron  interaction  physics  under  consideration,  we  turn  our  attention 
to  the  corresponding  transport  equation.  The  linear  Boltzmann  equation  for  the  angular  flux  of 
electrons  xp(r,  E,  H)  with  position  r(x,y,  z ),  energy  E,  and  direction  H  is  expressed  by 

i//(r,  E,  D)  =  Hb E,  D),  (10) 
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where  the  Boltzmann  collision  operator  is  defined  as 


n  00  n 

HBxp(r,E0 )  =  |  dE'  |  dtr[i;ei(r,E'  E0'  -h) 

Jo  J^n 

+  Zin(r,E'  E,n'  •■«)]  xp(r,E',0^  -  Zt(r,E)xp(r,E0), 


(ID 


and  the  total  interaction  cross  section  is  the  sum  of  the  elastic  and  inelastic  total  cross-sections, 
or  Zt(j,  E )  =  Zeii0(r,  E)  +  Zin  0(r,  E),  which  are  denoted  by  the  subscript  zero.  Eq.  (10)  is 
subject  to  the  following  boundary  condition 


y/(r,  E,  (p,/i  )  =xpb  (r,  E,  cp,  p  ),  r  E  dV,  fl  ■  n  <  0,  0  <  E  <  oo. 


(12) 


For  electrons,  elastic  and  inelastic  scattering  are  treated  separately  as  indicated  in  Eq.  (11). 
Furthermore,  it  is  assumed  that  elastic  scattering  occurs  without  energy  loss  and  angular 
deflection  from  inelastic  scattering  is  given  by  kinematics,  so  the  collision  operator  can  be 
expressed  as 

HBxp(r,E,n)  =  H^ir.E.D.)  +  H^f.E.O),  (13) 


where 


Hf^(r,E,n)=[  dn'Zel(r,EM' ■n)ip(r,E,tt)-Zel0(r,E)ip(r,E,n),  (14) 

-^47T 


and 


p  OO 

HBnxp(r,E0)  =  dE'  Zin(r,E'  E)  60'  -  f(E' ,QM(r,E' ,U) 
Jo 

-Zin, o(f,  E)xp{f,  E0). 


(15) 


Here  the  angle  of  the  primary  is 


Q  E  +  2  m0c2 
E  —  Q  +  2 m0c2  ’ 


(16) 


and  the  angle  of  the  secondary  is 


Q  E  +  2  m0c2 
E  —  Q  +  2 m0c2  ' 


(17) 


In  the  event  that  secondary  production  is  neglected,  deflection  of  the  primary  from  inelastic 
scattering  is  also  neglected  and  Eq.  (15)  reduces  to 
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(18) 


^  00 

H?nxj)(f,E,a)=  dE'  Zin(r,E'  ->  E)  \p{f  ,E'  —  Ein0(f ,  E)x/j(f ,  E  ,Tl). 

•'o 

As  previously  mentioned,  direct  numerical  solution  of  Eq.  (10)  is  not  realistic  for  most 
applications  of  interest  because  of  the  large  total  cross  sections,  EetQ  and  Ein0,  and  DCSs,  Eei 
and  Etn,  that  are  highly  peaked  about  small  deflections  and  energy  losses.  Because  the 
computational  effort  is  related  to  the  number  of  collisions  simulated  per  particle  history  and  the 
outcome  of  the  collision,  simulation  of  Coulomb  collisions  are  far  more  computationally 
intensive  than  the  interactions  characterizing  neutral  particles. 

4.2.  Moment-Preserving  Method 

In  the  MP  method,  a  reduced  order  physics  (ROP)  transport  equation  is  formed  by  replacing  the 
analog  DCSs  in  Eqs.  (14)  and  (18)  with  ROP  DCSs,  .  That  is, 

H*lip(f,E,£i)=  f  dn'Zel(r,E,n' -njip^E,^  -Zel/0(r,E)ip(r,E,n),  (19) 

■' 4n 

and 

^  CO 

H°nxp(r,E,n)  =  J  dE'  Ein(r,  E'  ->  E)  E\  fl)  -  Sinfi(f,  E)ip(r,  E,  fl).  (2Q) 

Although  simply  replacing  the  analog  DCS  with  an  ROP  DCS  may  seem  trivial  or  even  arbitrary, 
there  is  no  absence  of  rigor  in  this  method.  Particularly,  much  consideration  is  given  to  the  form 
and  properties  of  the  ROP  DCSs.  The  ROP  DCSs  are  constructed  such  that  they  are  smoother  or 
less-peaked  functions  of  deflection  angle  and  energy  loss  and  have  significantly  longer  mfps  than 
the  analog  DCSs.  Thus,  the  ROP  collision  operators  in  Eqs.  (11)  and  (12)  have  better  properties 
than  the  analog  collision  operators,  especially,  from  an  efficiency  standpoint. 

Beyond  efficiency,  there  are  additional  properties  of  the  ROP  collision  operators  that  set  this 
method  apart  from  other  approximate  methods.  For  example,  one  of  the  unique  characteristics  of 
this  method  is  that  the  integral  form  of  the  Boltzmann  collision  operators  are  maintained. 
Therefore,  the  description  of  the  underlying  transport  mechanics  is  not  lost,  specifically,  the 
correct  Markovian  feature  of  exponentially  distributed  collision  sites  [58].  Therefore,  special 
algorithms  for  handling  material  and  vacuum  interfaces  are  not  required.  Moreover,  exact 
treatment  of  collisions  as  Markov  processes  and  less-peaked  DCSs  with  longer  mfps  make  it 
practical  to  simulate  transport  with  a  single-event  method.  Implementation  of  single-event 
methods  is  very  straightforward  compared  to  other  methods  like  CH  that  is  considerably  more 
complicated.  In  fact,  Monte  Carlo  codes  with  pre-existing  single-event  algorithms  do  not  require 
any  retrofitting  when  implementing  the  MP  method,  because  this  method  treats  electrons  like 
neutral  particles. 


Approved  for  public  release;  distribution  is  unlimited. 

13 


The  MP  method  is  efficient  and  implementation  is  straightforward,  but  there  has  been  no 
mention  of  accuracy.  This  method  must  not  only  be  competitive  with,  and  potentially  superior  to 
CH  with  regard  to  efficiency  and  simplicity,  but  accuracy  as  well,  and  in  many  cases  it  is.  Unlike 
CH,  which  introduces  inherent  and  irreducible  limits  on  accuracy  as  a  result  of  the  underlying 
theory,  accuracy  is  systematically  controllable.  This  is  largely  a  result  of  the  moment-preserving 
strategies  that  are  central  to  this  approach.  The  moment-preserving  strategy  is  motivated  by 
Lewis  theory  [56,  57],  where  Lewis  proved  that  one  can  relate  space-angle  moments  of  the 
angular  flux  to  momentum-transfer  moments  of  the  elastic  scattering  DCS.  In  addition,  the 
eigenvalues  of  the  elastic  collision  operator  are  directly  dependent  on  the  momentum-transfer 
moments.  For  these  reasons,  it  is  prudent  to  construct  an  ROP  DCS  that  preserves  moments  of 
the  analog  DCS. 

Given  the  relationship  between  the  ROP  DCS  and  the  analog  DCS  moments,  the  following 
moment-preservation  constraints  are  a  natural  choice  when  constructing  an  ROP  DCS.  If  the 
analog  elastic  scattering  moments  are  given  by 

Zei,i(E)  =  2n  J  £1) 


the  moment  preserving  constraint  is 

Eei,i(E)  —  £ei,i(E),  l  —  1,2 

and  the  higher  order  moments  are  functions  of  the  lower  order  moments 

Zel,l(E)  =  f  (Z el, 1  Z el, 2>  ■■■’£ el, l)’ l  >  L. 


For  inelastic  scattering  the  moments  are  given  by 


rQmax 

ZinAE )  =  dQ  QjZin(E,Q), 


(22) 


(23) 


(24) 


and  the  moment  preserving  constraint  is  similar  and  given  by 

Zin,j (£)  =  Zin,j(E),j  —  1,2, 

Again,  the  higher  order  moments  are  functions  of  the  lower  order  moments 


(25) 


Ein,j(E )  —  Ein,2,  ■■■  ’  ZinA'j  >  /■ 


(26) 


By  constructing  an  ROP  DCS  that  preserves  moments  of  the  analog  DCS,  one  can  systematically 
control  the  accuracy  of  the  ROP  DCS  models.  That  is,  improvements  in  accuracy  are  achieved 
by  simply  preserving  more  moments  of  the  analog  DCS.  In  addition,  the  higher  order  moments 
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are  functions  of  the  exact  lower  order  moments,  so  the  higher  order  moments  are  good 
approximations  rather  than  being  neglected,  as  is  the  case  for  other  ROP  models  like  Fokker- 
Planck.  The  following  sections  present  the  two  ROP  DCS  models  that  are  used  in  this  paper. 
Given  the  discrete  ROP  DCS  model,  we  present  a  corresponding  derivation  of  the  ROP  transport 
model  for  a  discrete  elastic  and  inelastic  ROP  DCS  to  emphasize  the  difference  in  the  ROP 
transport  model  and  the  analog  model.  Lastly,  the  details  of  the  ROP  DCS  construction  process 
for  each  model  are  presented. 


4.2.1.  Reduced  Order  Physics  Differential  Cross-Sections 


In  this  section,  we  present  two  forms  of  the  elastic  and  inelastic  ROP  DCSs  that  are 
demonstrated  herein:  the  discrete  DCS  and  the  hybrid  DCS.  The  discrete  DCS  is  a  superposition 
of  discrete  points  and  weights.  One  of  the  benefits  of  the  discrete  DCS  is  the  simple  form  of  the 
DCS.  The  discrete  DCS  is  simple  to  sample  and  requires  significantly  less  memory  requirements 
than  DCS  data  because  only  a  few  points  and  weights  are  required  for  most  problems  of  inter-est. 
The  accuracy  and  efficiency  of  the  discrete  DCS  are  especially  promising  when  calculating 
integral  quantities  like  dose  [58].  We  define  the  discrete  DCS  for  elastic  scattering  as 


N 

Zel(?>E,n o)  =  <%0  -  <n]» 

n= 1 


and  for  inelastic  scattering  as 


N 

Zin(r,E,Q)  =  YjPn{E)8[Q  -Ynl 

n=l 


(27) 


(28) 


The  one  drawback  of  the  discrete  DCS  is  the  presence  of  discrete  artifacts  [61,  62],  especially,  if 
the  discrete  DCS  is  used  when  calculating  differential  quantities  in  thin  slabs.  However,  discrete 
artifacts  can  be  mitigated  by  use  of  the  hybrid  DCS,  while  still  achieving  efficiency  gains.  The 
hybrid  DCS  is  a  superposition  of  both  discrete  points  and  weights  and  a  smooth  function 
represented  by  an  analog  DCS.  In  previous  work  [61],  the  smooth  component  was  represented  by 
the  SR  DCS  over  [-1,  1].  The  screening  parameter  was  artificially  selected  such  that  the  smooth 
component  was  less  peaked  near  one.  Moments  of  the  smooth  component  are  then  subtracted 
from  the  analog  DCS  moments  and  this  difference  is  then  used  to  generate  the  discrete  scattering 
angles.  In  this  work,  a  slightly  different  representation  was  chosen  where  the  tail  is  represented 
exactly  by  the  analog  model  up  to  some  cutoff  point,  fi*0 .  Beyond  the  cutoff  point,  or  for  /r0  G 
[/To,  1],  a  discrete  representation  is  used.  The  resulting  hybrid  DCS  is 

N 

Eel(j,E,Ho)  =  Ze(E,Eo)  +  ~  ^ '  (29) 

n= 1 
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where  rf  (£,  /r0)  is  an  analog  DCS  for  /r0  G  [—1,  fi*0)  and  otherwise  zero.  The  cutoff,  fi*0,  is 
typically  chosen  to  be  near  the  peak  or  unity  to  gain  the  benefit  of  the  properties  of  the  discrete 
DCS,  while  capturing  the  large-angle  scattering  exactly  by  the  analog  DCS.  For  inelastic 
scattering,  the  cutoff,  Q,  is  selected  near  the  peak,  or  near  Qmin,  for  the  aforementioned  reasons. 

From  an  implementation  standpoint,  there  is  little  difference  between  the  discrete  and  hybrid 
DCS.  It  requires  the  ability  to  sample  hard  collisions  from  the  analog  DCS  and  soft  collisions 
from  the  discrete  DCS.  The  only  difference  in  generating  the  discrete  DCS  versus  the  hybrid 
DCS  is  that  the  moments  are  now  defined  over  a  partial  interval  corresponding  to  the  peak. 

Given  the  fonn  of  the  ROP  DCSs,  a  derivation  of  the  ROP  collision  operators  that  comprise  the 
ROP  transport  equation  is  presented. 


4.2.2.  Derivation  of  the  Reduced  Order  Physics  Collision  Operators 

The  ROP  DCS  is  constructed  such  that  the  singular  contribution  to  inscatter  and  outscatter  cancel 
(similar  to  the  FP  operator  [45]).  Ultimately,  the  purpose  of  constructing  such  a  DCS  is  that  the 
resulting  ROP  transport  equation  can  be  solved  accurately  and  efficiently  using  single-scatter 
models.  To  be  clear,  a  derivation  of  the  elastic  and  inelastic  ROP  collision  operators  is  presented. 
A  derivation  of  both  operators  is  presented  because  there  is  a  subtle  difference  between  the  two 
that  deserves  some  attention.  For  the  sake  of  simplicity,  the  discrete  DCS  is  used,  but  the  same 
ideas  carry  over  to  any  ROP  DCS. 


The  starting  point  is  the  elastic  collision  operator.  Substitution  of  Eq.  (27)  into  Eq.  (19)  gives 


N  +  l  riV  +  1 

HBif>(r,E,n)  =  J  ^an(E) 

4n  n= 1  Ln=l 


xp(f,  E,Tl).  (30) 


It  is  required  that  the  discrete  point  £w+1  =  1.  Now,  this  point  is  intentionally  separated  from  the 
remaining  N  points  and  weights  in  the  inscatter  and  outscatter  terms.  That  is, 

N 

HBelxp{r,E,n)  =  I  SnW-E.n')  - 

471  r  n=  1 

+  f  df2  aN+l(E)  g[Mo  _  (N]ip(f,E,Ur)  -  aN+1(E)ip(r,E,n). 

If  it  can  be  shown  that  the  last  two  tenns  in  Eq.  (31)  indeed  cancel,  the  resulting  ROP  transport 
equation  will  have  an  elastic  scattering  kernel  that  is  significantly  less  peaked  with  reduced  total 
cross  section  because  ccN+1  no  longer  contributes  to  the  total  cross  section.  It  is  not  difficult  to 
show  that  the  last  two  terms  cancel  because  /r0  =  fT  •  fl  =  1,  if  and  only  if  fT  =  fl  (see  ref.  [67] 
for  details).  Therefore,  the  following  is  true 


rv 

I 


an(E) 


ip(r,  E,U) 


no 
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f  an  aN+1<^E\[ii o  -  i]xp{r,E,H')  -  f  dn  Un+1^ s\a'  -a-  i]t/>(r,g,n') 

47T  ^7T  ^471 

=  ccN+1{E)xp(r,E,!l). 


(32) 


This  term  clearly  cancels  with  the  last  term  in  Eq.  (3 1).  However,  the  same  is  not  true  for  the 
inelastic  collision  operator,  which  is  now  shown. 

The  derivation  of  the  inelastic  ROP  collision  operator  is  easier  to  follow  after  a  change  of 
variables  change  variables  from  Q  to  ( E '  —  E).  This  gives 

N 

Sin(r,E,  Q)  =  ^  0n(E)  8[(Er  —  E)  —  Yn\.  (33 

n= 1 

Again,  Eq.  (33)  is  substituted  into  Eq.  (20)  and  the  argument  of  the  delta  function  is  rewritten  as 
E'  —  (E  +  Yn)  which  gives 


N  +  l 


H?nxl>(r,E,n)  =  J  dE'  ^  pn(E)  S[E'  —  (E  +  Yn)  ]  xf>(r,E',£l) 

0  n= 1 

rN+1  n 


n=  1 


4>(r,  E.D.). 


(34) 


Now,  the  singular  component  or  the  N+l  term  is  separated  from  the  inscatter  and  the  outscatter, 
the  inelastic  ROP  collision  operator  becomes 

_  r  a 

E,n)  =  J  dE'  2^  Pn(E)  S[E'  -(£  +  Yn)]  *(?,  E\  ft) 

0  n= 1 


^  Pn(E) 


Ln=l 

oo 


(35) 


4>(r,  E,  fl) 

/-  uu 

+  dE'pN+1(E)S[E'  ~(E  +  YN+1)]xl>(r,E',ci) 
3o 

-pN+1(E)ip(r,E,Cl). 


Carrying  out  the  integration  over  the  N+l  term  results  in  the  following  singular  contributions  to 
inscatter  and  outscatter 


pN+i(E)i/j(r,E  +  7W+I,n)  -  pN+1(E)ip(r,E,Q.). 


(36) 


The  terms  in  Eq.  (36)  do  not  cancel,  unless  Yn+1  —  0.  Typically,  the  singular  component  of  the 
inelastic  DCS  and  the  lower  bound  of  the  inscatter  or  energy-loss  moments  is  chosen  as  the  mean 
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ionization  potential.  Therefore,  yN+ 1  —  /Mey,  where  IMeV  is  the  mean  ionization  potential  in 
MeV.  So,  if  xp(r,  E  +  Yn+i>  is  expanded  as  a  Taylor  series  about  E,  one  can  get  a  sense  of  the 
error  introduced  by  choosing  Yn+i  —  ^Mev  •  That  is, 


pN+i(E)xp(r,E  +  yw+1,n)  -  (3N+1(E)\p{f,  E,D.) 


7  =  1  ' 


+  (3N+1(E)ip(f  ,E,Yl)  -  (3N+1(E)xl)(f ,  E ,Tl) 


(37) 


OO  j 

7  =  1  y' 

So,  to  first  order,  the  error  introduced  by  assuming  that  Yn+i  —  0  is  proportional  to  IMeV  which  is 
«  1 .  Therefore,  this  assumption  introduces  manageable  error. 


4.2.3.  Generation  of  the  Discrete  and  Hybrid  Differential  Cross-Sections 

The  procedure  for  constructing  both  the  discrete  and  hybrid  DCS  from  analog  DCS  moments  is 
described.  The  fonn  of  the  moment  preservation  constraints  in  Eqs.  (22)  and  (25)  is  unstable  to 
direct  numerical  inversion,  so  another  approach  similar  to  generation  of  Radau  quadrature  is 
taken  [68],  Initially,  attention  is  given  to  the  discrete  elastic  DCS  and  then  the  discussion  is 
extended  to  the  discrete  inelastic  DCS  and  the  hybrid  DCS. 

It  is  of  interest  to  obtain  a  DCS  that  satisfies  the  moment  constraint  in  Eq.  (22).  Given  the 
discrete  elastic  DCS  in  Eq.  (27),  a  system  of  equations  for  N  points  and  weights  (in  total  2N 
unknowns)  is  fonned.  Substitution  of  Eq.  (27)  into  the  right-hand-side  of  Eq.  (22)  results  in  the 
following  system  of  equations 

=  2tt  j 

N  1  1 

=  ^  «n(£)  j  -  <Tn]  (38) 

n= 1  -1 

N 

=  ^  tfnOOP^n) 

n=  1 

A  total  of  L  =  2N  equations  are  necessary  because  there  are  2N  unknowns.  That  is, 


Zei.i(E)  =  ai(£)f\0Ti)  +  ^(WiO^)  +  •"  + 

WE)  =  «i(£)W)  +  a2(E)P2«2)  +  -  +  a2N(E)P2((2N)  (39) 
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%el,2N(.E)  ~  (xi(.E')P2n(.^i)  +  al(£')^>2w(^2)  +  •"  +  a2w(^)^2iv(^2iv) 


The  system  fonned  in  Eq.  (39)  emphasizes  the  requirement  that  an  and  (n  are  obtained  such  that 
Legendre  moments  of  the  analog  DCS  are  preserved.  The  system  is  then  recast  into  one 
encountered  when  generating  Gauss-Radau  Quadrature  for  a  non-classical  weight  function  [69, 
70].  That  is, 


N 

Eel, l  =  ^  an(E)Pl(Zn)  +  «w+l(^)^((w+l  =  1)  -  (40) 

n= 1 

which  is  a  Gauss-Radau  Quadrature  system  for  a  non-classical  weight  function,  where  in  this 
case,  the  weight  function  is  the  analog  DCS.  Note  that  an  additional  unknown,  aN+1(E),  is 
added  in  Eq.  (40)  and  multiplied  by  Pi{^n+i  —  !)•  This  is  indicative  of  Radau  quadrature  and  an 
expression  for  detennining  aN+  l  (E)  is  given  below.  A  Radau  approach  is  selected  rather  than 
standard  Gauss  quadrature  because  Radau  ensures  that  one  point  will  correspond  to  the  peaked 
component  of  the  DCS  (that  is,  (N+1  =  1).  Once  the  discrete  points  and  weights  are  obtained,  the 
peaked  component  is  eliminated,  thus,  reducing  the  total  cross  section  after  renormalizing  the 
discrete  DCS.  This  is  equivalent  to  satisfying  the  moment  preservation  constraints  given  in  Eqs. 
(22)  and  (25). 


To  obtain  the  points  and  weights,  coefficients  of  monic  Legendre  polynomials  (a;  and  fij)  are 
mapped  to  the  coefficients  of  polynomials  orthogonal  to  the  analog  DCS  (a,  and  bj).  The 
algorithm  for  this  mapping  is  referred  to  as  the  modified  Chebyshev  algorithm  (MCA)  [70]  and 
requires  2N+1  moments  of  the  analog  DCS  and  2N+2  coefficients  of  monic  Legendre 
polynomials.  Given  a  successful  mapping  and  the  resulting  coefficients,  a,-  and  bj,  the  Golub  and 
Welsch  algorithm  [71]  is  used  to  obtain  the  eigenvalues  of  the  Jacobi  matrix.  The  Jacobi  matrix 
is  a  tridiagonal  matrix  where  the  diagonal  is  set  to  a,  and  the  off-diagonals  are  set  to  Vbj.  The 
eigenvalues  of  the  Jacobi  matrix  are  the  points  and  the  first  entry  of  each  corresponding 

eigenvector  squared  are  the  weights.  That  is,  (n  —  An(J)  and  an  —  (l/nl)  ,  where  Lisa 
eigenvector  matrix.  The  application  of  the  Golub  and  Welsch  algorithm  to  the  aforementioned 
Jacobi  matrix  will  result  in  Gauss  Quadrature  and  must  be  modified  according  to  Golub  [69]  for 
Radau  quadrature.  Therefore,  the  Jacobi  matrix  is  modified  such  that 


where 


Jn+i 


Jn 

bN  eTN 


b N  e. N 

aN+ 1  . 


—  i  u  Vn- i(1) 

“w+1  w^oT 


(41) 


(42) 
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Application  of  the  Golub  and  Welsch  algorithm  to  Eq.  (41)  will  result  in  N  +  1  points  and  N  +1 
weights  normalized  to  unity.  To  obtain  the  final  discrete  DCS,  the  N  +  1  point  and  weight  is 
eliminated  and  the  remaining  weights  are  then  scaled  by  the  analog  total  cross  section  or  Zei  0. 
The  total  cross  section  for  the  discrete  DCS  is  then 


■T el,0  —  ^ 


n=l 


(43) 


which  does  not  include  the  N  +  1  weight.  The  total  cross  section  in  Eq.  (43)  is  significantly 
reduced  depending  on  the  order  of  the  discrete  DCS,  the  particle  energy,  and  the  target  material, 
thus,  extending  the  mfp.  This  completes  the  process  of  generating  a  discrete  elastic  DCS. 


The  process  of  generating  a  discrete  inelastic  DCS  is  similar.  To  use  the  same  quadrature  tools, 
the  inelastic  DCS  must  be  mapped  to  an  elastic  DCS  because  the  bounds  on  the  elastic  DCS  are 
ideal  for  these  tools  (that  is,  [-1,  1]).  Given  a  mapping,  the  moments  of  the  inelastic  DCS  are 
related  to  Legendre  moments  of  an  ROP  elastic  DCS.  Points  on  (-1,  1]  are  generated  with 
corresponding  weights  and  then  mapped  back  to  (0,  Qmax ]  ■  The  mapping  is 

nr  \  Qmax  *. 

<200  =— 2~  (1-H).  (4, 


and  the  resulting  relationship  between  the  moments  is 


where 


i= 0 


(45) 


(46) 


This  summarizes  the  process  of  generating  the  discrete  elastic  and  inelastic  DCS.  Many  of  the 
same  ideas  carry  over  to  generation  of  the  discrete  component  of  the  hybrid  DCS. 


To  generate  the  discrete  points  and  weights  for  the  hybrid  DCS  a  cutoff  value  is  selected.  It 
should  be  selected  such  that  additional  accuracy  is  gained  while  still  maintaining  efficiency.  That 
said,  selection  of  the  cutoff  is  problem  dependent  and  mostly  a  heuristic  exercise. 


Given  a  cutoff,  the  following  moments  are  used  to  generate  the  discrete  points  and  weights 
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Zei.i  =  2n  I  dnPi(ji)ZelQi) 
■Vo 


(47) 


and 


’P„j=[  dQQ!Ein{Q). 

J  Q  min 


(48) 


In  both  cases,  to  use  the  DCS  generation  tools,  the  moments  must  be  mapped  to  the  appropriate 
domain  [-1,  1]  just  as  for  the  discrete  inelastic  DCS.  For  the  inelastic  hybrid  DCS  the  mapping 
does  not  change  significantly  from  Eqs.  (44)  and  (45)  and  is 


Q GO  =y  (!-  H), 

where  p  is  on  [-1,  1]  and  Q  is  on  [0,  Qcut].  Given  this  mapping,  the  moments  are  related  by 

*el’1  =  Z  C/'— Z  0  kz^k’ 


(49) 


j= o 


k= 0 


where  cj  is  given  by  Eq.  (46). 

The  map  for  the  hybrid  elastic  DCS  is  given  by 


(50) 


m'GO  =  — (l-iO  +  i, 


(51) 


where  is  on  p  [-1,  1]  and  ju'  is  on  [/i*0, 1],  Given  this  mapping,  the  moments  are  related  by 


Jel 


i=Yci 


,(-i  y 


j=o 


k= 0 


yD 

k  z' el.k ’ 


(52) 


where  b]k 


i  k 

bk  ~  j  ^/c(m)(1  -  ^  < 

_1  m= 0 


^ _ -Qm  2J+m+1 

m\  j  +  m  +  1  ’ 


(53) 


and  the  coefficient  is  given  in  Eq.  (46). 
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4.3.  Development  of  the  Geant4  Toolkit  Moment-Preserving  Method  Classes 

Reduced  order  physics  DCS  models  are  central  to  the  Moment-Preserving  method  and  therefore, 
the  majority  of  the  development  e  ort  is  in  generating,  storing,  and  accessing  the  ROP  DCSs. 
Geant4  considers  all  physical  interactions  as  processes,  requiring  implementation  of  the 
G4VEmProcess  class,  and  the  details  of  each  process  is  captured  by  the  model.  This  is 
accomplished  through  implementation  of  the  G4VEmModel  class,  which  is  the  primary  source 
of  development  for  the  MP  method  (assuming  the  ROP  DCS  are  available  in  the  form  of  a  DCS 
data  library).  The  reason  that  the  majority  of  the  development  occurs  when  implementing  the 
G4VEmModel  class  is  because  it  is  at  this  level  where  the  ROP  DCSs  are  stored  in  Geant4  data 
classes,  and  also  it  is  at  this  level  where  the  ROP  DCSs  are  accessed. 

It  should  noted  that  if  one  were  to  use  the  total  number  of  lines  of  code  required  to  implement  an 
electron  transport  algorithm  as  a  metric  for  quantifying  development  and  maintenance  effort,  the 
effort  associated  with  MP  method  is  about  ten  times  less  than  the  effort  required  by  the  Geant4 
default  multiple-scattering  model  used  to  transport  electrons.  That  is,  the  G4UrbanMscModel 
requires  roughly  2000  lines  of  code,  while  the  comparable  MP  elastic  model, 
G4DiscreteElasticModel,  only  requires  300  lines  of  code. 

The  remaining  sections  provide  detail  on  the  associated  physics  process  and  model  classes,  the 
cross-section  construction  classes,  and  the  cross-section  library  and  data  processing  tools. 

4.3.1.  Physics  Processes 

Physics  processes  for  discrete  and  hybrid  DCSs  were  implemented  and  include: 
G4DiscreteElasticProcess,  G4HybridSoftElasticProcess,  G4HybridHardElasticProcess, 
G4DiscreteInelasticProcess,  G4HybridSoftInelasticProcess,  and  G4HybridHardInelasticProcess. 
The  physics  process  is  relatively  simple  because  use  is  made  of  many  of  the  virtual  methods  in 
G4VEmProcess.  The  only  methods  that  required  implementation  include  a  method  to  initialize 
the  process  class  (this  consist  of  constructing  the  associated  physics  model  and  setting  a  few  data 
members)  and  a  method  that  determines  the  applicable  particles. 

4.3.2.  Physics  Models 

Physics  models  for  discrete  and  hybrid  DCSs  were  implemented  and  include: 
G4DiscreteElasticModel,  G4HybridSoftElasticModel,  G4HybridHardElasticModel, 
G4DiscreteInelasticModel,  G4HybridSoftInelasticModel,  and  G4HybridHardInelasticModel. 
Three  important  methods  are  included  in  the  physics  model  class:  an  initialization  method  where 
the  cross-section  data  is  read  in;  a  method  for  obtaining  the  total  cross-section;  and  a  method  for 
sampling  the  DCS.  Use  of  the  ROP  DCS  data  reduces  runtime  and  pre-existing  Geant4  data  tools 
were  utilized  that  reduced  development  overhead  and  eased  integration  of  the  models  into  the 
toolkit.  The  code  required  to  utilize  the  Geant4  data  classes  is  left  to  section  6.3,  which  describes 
the  use  the  Geant4  data  classes  in  greater  detail. 
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There  are  two  methods  remaining  for  the  physics  model  discussion,  a  method  for  obtaining  the 
ROP  total  cross-section  and  a  method  for  sampling  the  ROP  DCS.  Obtaining  the  total  cross- 
section  is  trivial  and  requires  a  table  look-up  as  a  function  of  energy  (some  minor  modifications 
to  the  Geant4  data  classes  were  required).  The  method  for  sampling  the  DCS  requires  two  steps. 
First,  the  energy  index  of  the  data  is  determined  using  the  previously  discussed  binary  search; 
however,  the  linear  interpolation  is  accomplished  through  random  sampling  using  the  weights 
from  a  linear  interpolation  on  log-log  scale  [64].  Given  an  ROP  DCS  evaluated  at  some  energy 
gird  point,  the  scattering  angle  is  detennined  through  inverting  a  discrete  CDF,  which  is  a  very 
simple  process. 

At  this  point,  no  mention  was  made  of  the  hybrid  cross-section.  As  seen  in  the  list  of  models,  the 
hybrid  cross-section  is  composed  of  two  models:  one  for  soft  collisions  and  one  for  hard 
collisions.  The  soft  collisions  are  given  by  the  discrete  cross  section  and  all  of  the  previous 
discussion  carries  over.  The  hard  collisions  are  given  by  the  analog  DCS.  For  the  partial-wave 
elastic  DCS  much  of  the  previous  discussion  also  carries  over.  However,  for  the  Moller  inelastic 
DCS,  an  analytical  expression  is  used  to  obtain  the  total  cross-section  and  a  rejection  technique  is 
used  to  sample  the  DCS. 

4.3.3.  Cross-Section  Library  and  Data  Processing 

A  ROP  cross-section  library  was  generated  for  the  partial-wave  and  Moller  DCS  for  1,  2,  4,  and 
8  points  and  weights.  The  libraries  are  fonnatted  such  that  the  Geant4  data  classes  could  be  used. 
For  each  material  and  number  of  points  and  weights,  there  are  two  data  les:  one  for  the  total 
cross-section  and  one  for  the  CDF.  The  les  are  named  accordingly.  For  example,  for  a  2-angle 
discrete  DCS  based  on  the  partial-wave  DCS  for  aluminum  the  two  les  are  named 
gbfp_pwe_tcs_13_2.dat  and  gbfp_pwe_cdf_13_2.dat,  where  the  first  number  is  the  atomic 
number  and  the  second  number  is  the  number  of  points  and  weights. 

There  are  two  primary  Geant4  data  classes  utilized  in  processing  and  storing  the  ROP  DCS  data: 
G4ElementData  and  G4PhysicsVector.  G4ElementData  is  a  very  powerful  class  that  stores  data 
for  a  particular  element  and  then  only  requires  the  atomic  number  to  retrieve  the  data.  Upon 
construction,  the  G4ElementData  object  requires  the  atomic  number  and  the  data  in  the  fonn  of  a 
pointer  to  a  G4Physics Vector  object,  so  the  G4ElementData  class  is  really  a  container  of 
G4PhysicsVector  objects.  The  data  is  stored  in  the  G4Physics Vector  object,  which  is  also  a  very 
powerful  class  because  as  long  as  the  data  is  properly  fonnatted,  one  must  simply  pass  the 
G4PhysicsVector  object  a  stream  of  the  data. 

Given  a  brief  description  of  the  implementation  of  the  MP  method  within  the  Geant4  toolkit,  we 
now  present  results. 
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5.  RESULTS  AND  DISCUSSION 


In  this  section,  a  wide  array  of  results  that  capture  the  key  features  of  the  Moment- 
Preserving  (MP)  method  is  presented.  In  particular,  the  key  features  of  this  method  demonstrated 
are:  systematic  accuracy,  efficient,  mathematically  robust,  versatile  and  simple. 

The  results  section  begins  by  demonstrating  the  first  feature  of  the  MP  method  through 
calculation  of  highly  differential  quantities  like  angular  distributions  and  energy  spectra.  In  these 
calculations,  the  MP  method  is  tested  under  the  strictest  possible  conditions  (that  is,  high-energy 
mono-energetic  pencil  beams  normally  incident  on  thin  slabs).  Under  these  conditions,  analog  or 
single-scatter  models  are  typically  required.  However,  it  is  shown  that  both  transmitted  and 
reflected  angular  distributions  and  energy-loss  spectra  can  be  resolved  through  use  of  the  hybrid 
DCS.  Though  the  emphasis  of  this  section  is  the  demonstration  of  the  systematic  nature  of  the 
MP  method,  efficiency  gains  of  at  least  five  times  analog  efficiencies  were  demonstrated  while 
maintaining  analog  level  accuracy  in  very  thin  slabs.  Under  these  extreme  conditions  where  the 
hybrid  model  is  successful,  the  discrete  model  tends  to  result  in  artifacts.  However,  it  is  possible 
to  utilize  the  discrete  model  in  thicker  slabs  where  the  benefit  of  efficiency  gains  is  significantly 
improved.  Following  the  angular  distribution  and  energy  spectrum  results,  longitudinal  and 
lateral  distributions  are  presented  where  the  first  two  features  of  the  method  are  again 
demonstrated  on  this  different,  but  important  quantity. 

Given  a  clear  understanding  of  the  systematic  feature  of  this  method,  results  for  less  extreme 
problem  conditions  are  presented  to  show  that  for  more  practical  applications  the  MP  method  is 
not  just  accurate,  but  also  very  efficient.  That  is,  we  show  that  for  1-D  and  2-D  dose  calculations 
the  MP  method  achieves  analog  level  accuracy  while  improving  efficiency  up  to  three  orders  of 
magnitude  over  analog  efficiencies.  We  show  that  this  is  true  for  low-Z  and  high-Z  materials,  for 
molecules  like  water  or  bone,  and  multi-region  problems.  Furthennore,  it  is  shown  that  material 
interfaces  in  multi-region  problems  do  not  introduce  additional  error  at  interfaces  as  does  the 
condensed  history  method.  This  is  because  the  MP  method  is  a  transport-  based  approximation 
and  the  benefit  of  this  type  of  approach  is  that  no  additional  algorithm  is  required  to  handle 
material  interfaces. 

Prior  to  this  work,  no  effort  had  been  made  to  validate  the  MP  method  through  comparison  with 
experimental  benchmarks.  Therefore,  several  results  are  presented  in  efforts  to  begin  this 
validation  process.  Specifically,  results  from  the  MP  method  are  compared  to  the  experimentally 
detennined  energy  deposition  profiles  (that  is,  the  Lockwood  data  [72]).  Similar  calculations  are 
made  with  the  Geant4  default  electromagnetic  physics  option  3,  so  that  accuracy  and  efficiencies 
for  the  MP  method  can  be  compared  to  the  Geant4  physics.  In  addition,  we  compare  charge 
deposition  results  generated  using  MP  method  with  experimentally  detennined  charge 
depositions  (that  is,  the  Tabata  data  [73]).  A  key  concept  to  point  out  in  comparing  with 
experimental  results  is  that  if  an  analog  model  exists  that  is  in  acceptable  agreement  with  an 
experimental  benchmark  and  moments  of  this  analog  model  are  readily  available,  one  can 
generate  reduced  order  physics  DCSs  based  on  the  aforementioned  analog  model  and  show 
similar  levels  of  agreement  while  significantly  improving  efficient.  In  this  sense,  the  versatility 
and  simplicity  of  the  method  is  demonstrated.  That  is,  to  improve  agreement  through  use  of  a 
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different  analog  model  algorithmic  changes  are  not  required;  one  must  simply  obtain  moments  of 
some  preferable  analog  model,  generate  an  ROP  DCS  library  for  this  analog  model,  and  run  the 
calculation. 

Finally,  for  completeness  the  MP  method  is  applied  to  a  space  weather  application 
to  show  that  this  method  is  not  just  effective  for  theoretical  calculations.  In  particular,  the  total 
response  function  is  generated  for  the  CEASE  detector  telescope  using  an  analog  model,  a 
discrete  model,  and  the  default  Geant4  electromagnetic  physics  with  option  3. 

In  all  comparisons,  the  analog  benchmark  is  obtained  by  using  an  analog  elastic 
DCS  and  an  analog  inelastic  DCS  (with  exception  of  the  validation  section).  The  analog 
benchmark  is  numerical  rather  than  experimental,  so  to  some  degree  it  is  idealized.  This  type  of 
benchmarking  is  required  to  illustrate  how  accuracy  is  achieved  through  preservation  of  the 
analog  DCS  moments.  In  the  following  sections,  accuracy  and  efficiency  is  measured  with 
respect  to  the  analog  benchmark. 

5.1.  Angular  Distributions  and  Energy  Spectra 

In  this  section,  the  accuracy  of  the  MP  method  is  tested  under  conditions  that  are  often  times 
impractical  to  simulate  without  the  use  of  an  analog  model.  That  is,  the  transport  of  a  mono- 
energetic  pencil  beam  of  electrons  with  energies  above  several  hundred  keV  in  thin  slabs.  Of 
particular  interest  is  the  calculation  of  reflected  and  transmitted  angular  distributions  and  energy 
spectra  in  slabs  with  varying  thick-  nesses  down  to  100  analog  elastic  mfps  (~  1  to  100  pm).  It  is 
shown  that  under  these  conditions  the  MP  method  is  effective  at  resolving  angular  distributions 
and  energy  spectra  through  the  use  of  suitable  reduced  order  physics  (ROP)  DCSs.  Under  these 
extreme  conditions,  a  hybrid  DCS  is,  in  most  cases,  required  to  resolve  these  distributions.  While 
it  is  possible  to  resolve  highly  peaked  distributions  with  analog  level  accuracy  using  the  MP 
method,  this  is  typically  accompanied  by  losses  in  efficiency.  Nonetheless,  the  ability  to 
systematically  control  accuracy  such  that  one  can  predict  angular  distributions  and  energy 
spectra  for  highly-peaked  scattering  in  thin  slabs  is  a  strong  feature  of  this  method. 

While  it  is  true  that  one  cannot  expect  to  realize  significant  efficiency  gains  with  analog  level 
accuracy  under  the  aforementioned  conditions,  in  more  realistic  settings  (that  is,  thicker  regions) 
it  is  possible  to  relax  the  ROP  models  such  that  both  analog  level  accuracy  and  significant 
efficiency  gains  are  achieved.  Relaxation  of  the  ROP  models  is  possible  in  thicker  slabs  because 
the  initial  pencil  beam  experiences  more  spreading  in  space,  angle,  and  energy  in  thicker  slabs. 
This  is  simply  an  effect  of  the  number  of  collisions  sustained  by  an  electron  while  traversing  a 
medium.  In  thicker  slabs,  electrons  suffer  more  collisions;  thus,  causing  additional  spreading  of 
the  initial  state  of  the  beam.  With  additional  spreading  of  the  beam,  less  information  in  the  fonn 
of  analog  DCS  moments  is  required  to  resolve  angular  distributions  and  energy  spectra. 
Therefore,  the  ROP  DCS  can  be  relaxed  or  models  preserving  fewer  moments  can  be  utilized. 

In  the  following  two  sections,  the  impact  of  the  size  of  the  slab  on  the  accuracy  of  the  MP 
method  is  demonstrated.  Both  discrete  and  hybrid  models  are  tested  for  angular  distribution  and 
energy  spectrum  calculations  in  low-Z  and  high-Z  slabs  with  varying  thicknesses.  The  thickness 
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of  the  slab  is  measured  with  respect  to  the  analog  elastic  mfp  corresponding  to  the  source  particle 
energy  and  the  target  material.  Results  for  slabs  with  thicknesses  of  100,  300,  1000,  and  3000 
analog  elastic  mfps  are  presented  (in  the  remaining  discussion  mfp  implies  analog  elastic  mfp). 

In  Figure  3,  problem  setup  is  described. 


5.1.1.  Angular  Distributions 

Reflected  and  transmitted  angular  distributions  are  presented  below  for  one-dimensional  slabs 
composed  of  aluminum  or  gold  with  thicknesses  of  100,  300,  1000,  and  3000  mfps.  The  source  is 
positioned  at  x  =  0  with  a  direction  of  Q  =  (1,  0,  0).  A  total  of  4E+07  source  particles  are 
simulated  when  calculating  the  angular  distributions.  The  analog  benchmark  is  a  solution  to  the 
aforementioned  problem  using  analog  Monte  Carlo,  where  elastic  scattering  is  given  by  the 
partial- wave  DCS  and  inelastic  scattering  is  given  by  the  Moller  DCS.  Uncertainties  associated 
with  these  results  are  within  1%  in  most  bins,  such  that  one  can  state  conclusively  that  good 
agreement  exists  between  the  ROP  models  and  the  analog  benchmark. 

First,  the  most  challenging  problem  is  presented.  That  is,  calculation  of  transmitted  angular 
distributions  for  10000-keV  electrons  incident  on  an  aluminum  or  gold  slab  100  mfps  thick.  In 
Figure  4,  transmitted  angular  distributions  in  aluminum  and  gold  computed  using  the  discrete 
and  hybrid  DCSs  are  compared  to  the  analog  benchmark.  There  are  a  few  features  to  note  in 
Figure  4.  The  peakedness  of  this  distribution  is  extreme  and  varies  about  three  orders  of 
magnitude  over  only  10  degrees.  This  level  of  peakedness  is  difficult  to  resolve  with  a  discrete 
elastic  DCS  and  results  in  the  discrete  artifacts  seen  clearly  in  Figure  4.  The  difficulty  in 
resolving  highly  peaked  distributions  using  the  discrete  DCS  results  from  the  form  of  the  DCS. 
That  is,  electrons  can  only  scatter  through  N  discrete  angles  detennined  by  the  order  of  the  DCS. 
Therefore,  the  N  discrete  angles  are  favored  in  the  angular  distribution  because  in  thin  slabs 
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electrons  do  not  suffer  enough  collisions  such  that  various  combinations  of  scattering  events 
smooth  out  the  artifacts.  However,  through  use  of  the  hybrid  model  the  discrete  artifacts  are 
mitigated  and  the  only  noticeable  differences  in  the  hybrid  model  solution  and  the  analog 
benchmark  are  in  the  tail  where  the  differences  in  the  solutions  are  statistically  insignificant.  In 
Figure  4,  the  impact  of  target  atomic  number  is  shown,  where  for  increasing  Z  the  distribution  is 
less  peaked.  However,  the  impact  of  the  atomic  number  on  the  peakedness  of  the  scattering  is  not 
significant  enough  to  dramatically  improve  the  discrete  results  in  thin  slabs  for  10000-keV 
electrons. 


Figure  4:  Transmitted  angular  distributions  for  10000-keV  electrons  on  aluminum  and 
gold. 

Even  in  aluminum  slabs  with  thicknesses  of  3000  mfps  for  10000-keV  electrons,  the  discrete 
DCS  results  in  artifacts  as  seen  in  Figure  5.  This  is  an  indication  of  the  extreme  peakedness  of 
the  scattering  at  higher  energies  and  in  this  regime  a  hybrid  DCS  is  required  to  resolve  the 
transmitted  angular  distribution  unless  additional  angles  are  used.  However,  it  is  clear  by  Figure  5 
that  in  thicker  slabs  where  particles  undergo  thousands  of  collisions  that  the  discrete  artifacts  are 
greatly  reduced.  Furthermore,  the  impact  of  the  atomic  number  is  seen  in  Figure  5,  where  the 
discrete  artifacts  are  less  pronounced  in  the  gold  slab  because  scattering  of  electrons  by  high-Z 
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materials  is  less  peaked.  Nonetheless,  it  is  always  possible  to  utilize  a  hybrid  model  to  resolve 
angular  distributions  overwhelmed  by  discrete  artifacts. 

In  thicker  slabs,  it  is  possible  to  relax  the  cutoff  to  p*  =  0.99  for  the  hybrid  model  improving  the 
efficiency  of  the  calculation  while  remaining  accurate.  The  results  in  Figure  4  and  Figure  5 
indicate  that  it  is  possible  to  resolve  angular  distributions  in  highly-peaked  scattering  regimes  by 
systematically  increasing  the  accuracy  of  the  ROP  DCS  through  preservation  of  additional 
moments.  Of  course,  increasing  accuracy  will  reduce  the  efficiency  of  the  calculation,  but  under 
these  conditions  (highly-peaked  scattering  in  thin  slabs)  analog  Monte  Carlo  efficiencies  are 
typically  manageable,  so  efficiency  gains  of  two  to  five  times  faster  than  analog  Monte  Carlo  is 
considered  a  significant  improvement.  As  was  pointed  out,  the  emphasis  of  this  section  was  not 
to  necessarily  demonstrate  orders  of  magnitude  efficiency  gains,  but  rather  to  show  that  the  ROP 
models  limit  to  analog  level  accuracy  even  under  extreme  conditions.  That  said,  it  is  of  interest  to 
maximize  efficiency  gains  whenever  possible.  Therefore,  the  following  results  provide  a  sense  of 
the  accuracies  associated  with  a  less-extreme  scattering  regime.  Efficiency  results  are  presented 
later  in  section  5.1.3  and  indicate  that  it  is  possible  to  resolve  angular  distributions  efficiently  (up 
to  two  orders  of  magnitude  more  efficient  than  analog)  and  accurately,  especially,  in 
less-extreme  scattering  regimes. 


Figure  5:  Transmitted  angular  distributions  for  1000-keV  electrons  on  aluminum  and  gold. 


The  following  figures  present  transmitted  angular  distributions  in  a  less-peaked  scattering 
regime.  That  is,  1000-keV  electrons  incident  on  gold  slabs  with  varying  thicknesses.  As  noted, 
the  peakedness  of  the  scattering  is  a  function  of  particle  energy  and  the  target  atomic  number. 
With  decreased  particle  energy  and  increased  atomic  number,  the  peakedness  is  reduced. 
Nonetheless,  even  for  1000-keV  electrons  on  gold  the  problem  is  still  extremely  anisotropic  with 
respect  to  neutral  particle  scattering.  In  Figure  6,  the  impact  of  slab  thickness  and  in  turn,  the 
effectiveness  of  the  discrete  model  is  demonstrated.  In  Figure  6  discrete  artifacts  are  present  for 
slabs  100  and  300  mfps  thick,  but  the  hybrid  model  is  in  good  agreement  in  these  cases. 
However,  in  Figure  6,  discrete  models  with  at  most  4-angles  are  sufficient  when  resolving  the 
transmitted  angular  distribution.  In  fact,  in  Figure  6  the  discrete  artifacts  resulting  from  a  single¬ 
angle  discrete  model  are  almost  negligible  and  though  there  is  not  perfect  agreement  the  general 
behavior  of  the  transmitted  angular  distribution  is  captured. 
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Figure  6:  Impact  of  slab  thickness  on  transmitted  angular  distributions  on  1000-keV 
electrons  on  gold. 


We  now  present  reflected  angular  distributions  for  1000-keV  and  10000-keV  electrons  on 
aluminum  or  gold  slabs  with  thicknesses  of  100,  300,  1000,  and  3000  mfps.  In  Figure  7  reflected 
angular  distributions  for  10000-keV  electrons  on  gold  are  presented.  First,  note  the  distributions 
in  Figure  7  are  significantly  reduced  in  magnitude  relative  to  the  transmitted  angular 
distributions  for  10000-keV  electrons.  For  highly-  peaked  scattering,  a  very  small  fraction  of 
particles  are  reflected.  Slab  thickness  has  a  similar  impact  on  reflected  distributions  as  for 
transmitted  distributions.  That  is,  with  increasing  slab  thickness  electrons  suffer  more  collisions 
before  being  reflected,  spreading  the  distributions  in  angle.  In  general,  the  discrete  model  tends 
to  have  the  correct  behavior;  however,  the  distribution  is  roughly  two  to  five  times  greater  in 
magnitude.  Once  again,  the  hybrid  model  can  be  used  in  all  cases.  The  disagreement  between  the 
analog  benchmark  and  the  hybrid  models  is  statistical  because  very  few  particles  are  reflected. 
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Figure  7:  Impact  of  slab  thickness  on  reflected  angular  distributions  for  10000-keV 
electrons  on  aluminum. 

In  Figure  8  angular  distributions  for  1000-keV  electrons  incident  on  gold  slabs  with  thicknesses 
of  100,  300,  1000,  and  3000  mfps  are  presented.  For  lower  energies  in  high-Z  materials  the 
scattering  is  less  peaked  and  the  ROP  model  can  be  relaxed  under  these  conditions.  Specifically, 
reflected  angular  distributions  generated  using  the  discrete  model  are  not  overwhelmed  by 
artifacts  as  seen  in  Figure  8.  For  slabs  of  sufficient  thickness,  a  hybrid  model  is  not  required  and 
the  more  efficient  discrete  model  can  be  utilized.  For  slabs  that  are  3000  mfps  thick  (Figure  8),  a 
single-angle,  single-energy  discrete  model  provides  noteworthy  agreement. 
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Figure  8:  Impact  of  slab  thickness  on  reflected  angular  distributions  for  1000-keV  electrons 
on  gold. 

5.1.2.  Energy  Spectra 

Next,  reflected  and  transmitted  energy  spectra  are  examined.  The  simulation  characteristics, 
including  slab  thickness,  material  types,  and  the  number  of  source  particles,  are  the  same  as 
described  in  the  previous  section.  Discrete  inelastic  DCS  models  are  not  presented  because  a 
sufficient  number  of  inelastic  collisions  do  not  occur  in  thin  slabs  overwhelming  the  spectra  with 
discrete  artifacts.  Therefore,  the  focus  is  on  two  different  hybrid  inelastic  DCS  models.  The 
hybrid  models  include  DCSs  with  Q*  =  10-keV  and  one  or  two  discrete  energies.  In  each  of  the 
following  results,  elastic  scattering  is  modeled  by  a  discrete  four-angle  DCS.  Again,  the  most 
challenging  problem,  10000-keV  electrons  on  aluminum  slabs,  is  considered  first.  In  Figure  9, 
transmitted  energy-loss  spectra  are  presented  for  aluminum  slabs  with  thicknesses  of  100,  300, 
1000,  and  3000  mfps.  As  seen  in  Figure  9,  it  is  possible  to  resolve  the  transmitted  energy-loss 
spectra  with  a  sufficiently  accurate  hybrid  inelastic  model.  For  10000-keV  electrons  on 
aluminum,  a  two-energy  hybrid  model  is  required.  How-  ever,  for  1000-keV  electrons  on  gold,  it 
is  possible  to  relax  the  inelastic  model  to  a  single  energy  with  the  same  cutoff.  In  Figure  10, 
results  for  1000-keV  electrons  on  gold  are  presented.  Under  these  conditions,  inelastic  scattering 
is  also  less  peaked  and  the  hybrid  model  can  be  relaxed.  However,  a  discrete  representation  is 
still  not  sufficient. 
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Figure  9:  Impact  of  slab  thickness  on  transmitted  energy-loss  spectra  for 
electrons  on  aluminum. 


10000-keV 
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Figure  10:  Impact  of  slab  thickness  on  transmitted  energy-loss  spectra  for  1000-keV 
electrons  on  gold. 


Reflected  energy-loss  spectra  for  10000-keV  and  1000-keV  electrons  on  gold  are  presented  next. 
Reflected  energy-loss  spectra  for  10000-keV  electrons  on  aluminum  are  not  presented  because 
the  reflected  electrons  experience  almost  no  inelastic  collisions.  This  is  true  even  for  gold  as  seen 
in  Figure  11  for  slabs  that  are  100  and  300  analog  mfps  thick,  where  the  distributions  are  nearly 
singular  about  zero  energy-loss.  That  said,  the  hybrid  model  also  predicts  this  nearly  singular 
behavior.  For  thicker  slabs,  1000  and  3000  analog  mfps,  the  reflected  energy-loss  spectra  spreads 
out  more.  Flere,  the  hybrid  model  is  in  good  agreement  in  the  statistically  significant  regions  of 
the  spectra  (that  is,  the  peaked  region  of  the  spectra). 
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Figure  11:  Impact  of  slab  thickness  on  reflected  energy-loss  spectra  for  10000-keV 
electrons  on  aluminum. 
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In  Figure  12,  reflected  energy-loss  spectra  are  presented  for  1000-keV  electrons  on  gold.  Once 
again,  in  this  less-peaked  regime  the  spectra  are  not  nearly  as  singular  and  a  relaxed,  single¬ 
energy  hybrid  model  can  be  utilized. 


Figure  12:  Impact  of  slab  thickness  on  reflected  energy-loss  spectra  for  1000-keV  electrons 
on  gold. 

5.1.3.  Efficiencies  for  Thin  Slab  Problems 

In  this  section,  efficiency  results  for  the  previously  discussed  thin  slab  problems  are  presented  in 
Table  1  and  Table  2  for  1000-keV  and  10000-keV  on  aluminum  or  gold  slabs  of  various 
thicknesses  respectively.  In  general,  efficiency  gains  depend  on  the  slab  thickness,  the  target 
atomic  number,  and  the  energy  of  the  particle.  The  dependence  on  the  problem  geometry,  or  the 
slab  thickness  in  particular,  is  captured  in  the  following  tables.  The  greatest  efficiencies  are 
realized  in  thicker  slabs  where  particles  undergo  more  collisions.  In  thin  slabs  (100  and  300 
analog  mfps),  efficiency  gains  range  from  5  to  60  times  faster  than  analog  Monte  Carlo.  Whereas 
for  thicker  slabs,  efficiency  gains  range  from  one  to  three  orders  of  magnitude  faster  than  analog 
Monte  Carlo. 

In  previous  sections,  it  was  mentioned  that  there  is  a  trade-of  between  accuracy  and  efficiency. 
This  will  be  clarified  in  remainder  of  this  discussion.  First,  remember  that  under  the  most 
extreme  conditions  (10000-keV  electrons  on  100  to  300  mfps  of  aluminum)  hybrid  models  were 
required  to  resolve  angular  distributions  and  energy  spectra.  The  efficiency  gains  associated  with 
these  models,  under  these  conditions  range  from  15-47  times  the  efficiency  of  an  analog  Monte 
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Carlo  calculation.  What  would  take  a  day  now  requires  only  an  hour  without  sacrificing 
accuracy.  Now,  note  that  it  was  possible  to  resolve  angular  distributions  for  1000-keV  electrons 
on  gold  slabs  with  thicknesses  of  1000  and  3000  mfps  using  various  discrete  models.  Under 
these  conditions,  efficiency  gains  of  150-500  times  the  efficiency  of  an  analog  Monte  Carlo 
calculation  are  realized  without  sacrificing  accuracy.  In  this  case,  what  would  take  a  day  takes 
less  than  ten  minutes. 

Table  1:  Efficiency  gains  for  various  ROP  DCSs  when  simulating  1000-keV  and  10000-keV 
electrons  incident  on  aluminum  slabs  100,  300, 1000,  and  3000  mfps  thick. 


Reduced  Order  Physics  Model 

Slab  Width 

Particle 

1 -Angle 

4-Angles 

4-Angles 

1 -Angle 

Mo  =  0.99 

4-Angles 

1 -Energy 

(mfps) 

Energy 

1 -Energy 

1 -Energy 

4-Energies 

4-Energies 

Q*  =  10  keV 

100 

1000-keV 

21 

18 

15 

5 

10 

10000-kev 

22 

22 

21 

16 

13 

300 

1000-keV 

56 

43 

33 

16 

22 

10000-kev 

60 

62 

55 

48 

27 

1000 

1000-keV 

174 

63 

56 

11 

41 

10000-kev 

200 

170 

145 

86 

85 

3000 

1000-keV 

475 

178 

109 

21 

76 

10000-kev 

610 

556 

360 

114 

71 

Table  2:  Efficiency  gains  for  various  ROP  DCSs  when  simulating  1000-keV  and  10000-keV 
electrons  incident  on  gold  slabs  100,  300, 1000,  and  3000  mfps  thick. 


Reduced  Order  Physics  Model 

Slab  Width 

Particle 

1 -Angle 

4-Angles 

4-Angles 

1 -Angle 

Mo  =  0.99 

4-Angles 

1 -Energy 

(mfps) 

Energy 

1 -Energy 

1 -Energy 

4-Energies 

4-Energies 

Q*  =  10  keV 

100 

1000-keV 

20 

16 

15 

8 

12 

10000-kev 

19 

20 

19 

17 

20 

300 

1000-keV 

61 

37 

33 

10 

26 

10000-kev 

62 

57 

54 

86 

45 

1000 

1000-keV 

174 

63 

56 

18 

30 

10000-kev 

200 

170 

145 

89 

39 

3000 

1000-keV 

565 

134 

122 

23 

91 

10000-kev 

1174 

861 

591 

255 

254 

In  section  5.1.1  and  5.1.2,  it  was  shown  that  for  a  given  ROP  DCS  model,  the  accuracy  of  the 
result  depends  on  the  peakedness  of  the  distribution  being  resolved.  The  peakedness  of  the 
distribution,  in  turn,  depends  on  the  slab  thickness  and  the  regime  of  scattering,  which  is  a 
function  of  particle  energy  and  target  atomic  number.  Regardless,  an  ROP  DCS  can  be  made  to 
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preserve  additional  moments  of  the  analog  DCS;  thus,  refining  the  model  and  achieving  analog 
level  accuracy  under  very  extreme  simulation  conditions.  While  it  is  not  possible  to  achieve 
several  orders  of  magnitude  efficiency  gains  under  these  conditions,  the  ROP  models  were  at 
least  five  times  faster  than  analog  Monte  Carlo  and  up  to  45  times  faster  in  some  cases.  Under 
more  relaxed  conditions,  up  to  three  orders  of  magnitude  efficiency  gains  were  achieved. 
Regardless,  accuracy  and  efficiency  suitable  for  a  wide  variety  of  problems  can  be  realized  by 
adjusting  the  ROP  DCSs. 

5.2.  Longitudinal  and  Lateral  Distributions 

At  this  point,  it  is  clear  that  accuracy  and  efficiency  is  problem  dependent.  Nonetheless,  it  is 
possible  to  refine  the  ROP  DCS  such  that  sufficient  levels  of  accuracy  and  efficiency  are 
realized.  In  this  section,  a  few  additional  results  are  presented  that  overlap  with  the  previous 
section  in  the  sense  that  problems  corresponding  to  scattering  regimes  ranging  from  extreme 
peakedness  to  less-extreme  peakedness  are  tested.  Here,  results  in  connection  with  Lewis  theory 
are  presented  to  demonstrate  the  moment -preserving  property  of  this  method,  while  again 
demonstrating  the  accuracy  of  this  method  is  systematically  controllable. 

In  particular,  longitudinal  and  lateral  distributions  for  100-keV,  1000-keV,  and  10000-keV 
electrons  after  traveling  a  distance  of  100,  300,  1000,  and  3000  mfps  in  an  infinite  medium  of 
copper  are  presented  (similar  to  Benedito  et  al.  [64]).  In  these  problems,  energy-loss  is  not 
considered.  The  longitudinal  and  lateral  distributions  are  generated  using  the  analog  Monte  Carlo 
method  where  elastic  scattering  is  given  by  the  partial-wave  DCS.  These  distributions  are 
referred  to  as  the  analog  benchmark  and  are  compared  with  several  solutions  generated  using 
discrete  and  hybrid  models.  In  all  of  the  results  presented  in  this  section,  4E+07  electrons  were 
simulated  for  each  model  and  uncertainties  associated  with  the  majority  of  the  results  are  within 
1%;  especially,  in  the  highly  probable  regions.  However,  in  regions  where  the  distribution  is 
small  with  respect  to  the  maximum  value  (for  example,  in  the  tails  of  the  distributions),  the 
results  are  statistically  insignificant. 

In  Figure  13,  a  diagram  of  the  problem  simulated  is  presented.  The  electron  starts  at  s  =  0  with 
an  initial  direction  and  travels  until  reaching  a  distance  of  s  =  smax .  At  this  point,  the  electrons 
longitudinal  displacement,  or  the  projection  of  the  path  traveled  onto  the  initial  trajectory,  and 
lateral  displacement,  or  the  orthogonal  projection  of  the  pathlength,  is  tallied.  Although  it  is  not 
clear  from  Figure  13,  it  is  possible  for  an  electron  to  turn  around  and  travel  in  directions 
opposite  to  the  initial  direction.  In  these  cases,  it  is  possible  for  the  electron  to  have  a  negative 
longitudinal  displacement.  The  same  is  not  true  for  lateral  displacement  because  the  lateral 
displacement  is  a  measure  of  radius  and  therefore,  non-negative. 
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S=100 


(x2+y2)I/2 


In  addition  to  the  longitudinal  and  lateral  distributions,  a  few  Lewis  moments  are  compared 
including  (z)  and  (x2  +  y2)  or 


/i  /i  OO  /i  OO  /i  oo 

(z)  —  df]  dx  dy  dz  zxp(x,y,z,s,n) . 

J  47T  *•'  —  CO  J  —  OO  j  —  OO 


(54) 


and 


(x2  +  y2)  =  f 

J  4 


dn 


An 


r>  OO  OO  C 

I  dx  I  dy  I 

j  —  CO  J  —  CO  J  — 


dz  (x2  +  y2)t/>(x,  y,  z,  s,  fi)  . 


(55) 


Here,  the  Monte  Carlo  method  is  used  to  carry  out  the  integrals  in  Eqs.  (54)  and  (55),  by 
simulating  the  particle  transport  and  tallying  the  longitudinal  and  lateral  displacement  after 
traveling  100,  300,  1000,  or  3000  analog  mfps.  These  results  are  presented  in  the  following 
tables  for  each  energy.  As  predicted  by  Lewis  theory  and  seen  in 

Table  3,  models  preserving  at  least  Ze[1  will  preserve  (z).  Therefore,  even  the  very  efficient 
single-angle  model  will  have  the  correct  average  longitudinal  displacement.  Once  again  as 
predicted  by  Lewis  theory  and  seen  in  Table  4,  models  preserving  at  least  Eei  l  and  L'e^will 
preserve  (x2  +  y2).  Therefore,  the  single-angle  model  that  preserves  Zet  l  and  Zet  2will  have  the 
correct  average  lateral  displacement  as  well.  Models  preserving  additional  moments  are  not 
presented  in 

Table  3  and  Table  4  because  these  results  are  redundant.  Preservation  of  average  longitudinal 
and  lateral  displacement  are  important  to  electron  transport  methods  and  in  many  cases  these 
methods  seek  to  preserve  at  least  average  longitudinal  and  lateral  displacement.  One  of  the  major 
distinctions  between  condensed  history  and  the  MP  method  is  that  typically  in  condensed  history 
the  underlying  multiple  scattering  theory  is  only  setup  to  preserve  Zel  l  and  Zel  2,  while  the  MP 
method  can  preserve  an  arbitrary  number  of  Legendre  moments  guaranteeing  preservation  of 
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higher-order  Lewis  moments.  Not  to  mention,  the  most  simple,  efficient  ROP  model  associated 
with  the  MP  method  preserves  at  least  the  average  longitudinal  and  lateral  distances. 

Table  3:  Average  longitudinal  displacement  for  100-keV,  1000-keV,  and  10000-keV 
electrons  in  copper  traveling  a  distance  of  100,  300, 1000,  and  3000  mfps. 


Pathlength 

(mfps) 

Particle  Energy 
(keV) 

(z) 

Analog 

1 -Angle 

Rel.  Unc. 

100 

0.8432 

0.8433 

0.00006 

100 

1000 

0.986667 

0.98667 

0.00001 

10000 

0.999644 

0.999644 

0.000001 

100 

0.61815 

0.61820 

0.0002 

300 

1000 

0.96069 

0.96070 

0.00004 

10000 

0.998931 

0.998932 

0.000005 

100 

0.2761 

0.2763 

0.001 

1000 

1000 

0.87673 

0.87680 

0.0001 

10000 

0.996449 

0.996448 

0.00001 

100 

0.09493 

0.09497 

0.006 

3000 

1000 

0.6861 

0.6863 

0.0007 

10000 

0.989399 

0.989387 

0.00005 

Table  4:  Average  lateral  displacement  for  100-keV,  1000-keV,  and  10000-keV  electrons  in 
copper  traveling  a  distance  of  100,  300, 1000,  and  3000  mfps. 


Pathlength 

(mfps) 

Particle  Energy 
(keV) 

(x2 

Analog 

+  y2) 

1 -Angle 

Rel.  Unc. 

100 

0.154457 

0.154456 

0.0003 

100 

1000 

0.016420 

0.016415 

0.0008 

10000 

0.000455 

0.000454 

0.003 

100 

0.26942 

0.26944 

0.0004 

300 

1000 

0.04683 

0.04682 

0.0009 

10000 

0.001365 

0.001364 

0.004 

100 

0.2365 

0.2366 

0.0007 

1000 

1000 

0.1315 

0.1314 

0.001 

10000 

0.004524 

0.004525 

0.0009 

100 

0.11012 

0.11016 

0.002 

3000 

1000 

0.2532 

0.2531 

0.001 

10000 

0.01338 

0.01339 

0.004 

Although  preservation  of  Zei  l  and  2  guarantees  preservation  of  (z)  and  ( x 2  +  y2),  resolving 
the  longitudinal  and  lateral  distributions  requires  preservation  of  additional  moments.  Similar  to 
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results  from  the  previous  section,  discrete  artifacts  impact  the  shape  of  the  longitudinal  and 
lateral  distributions.  The  impact  of  the  artifacts  is  again  dependent  on  the  peakedness  of  the 
scattering  and  the  number  of  collisions  suffered  by  the  electrons  before  tallying  their 
displacement. 

In  Figure  14  through  Figure  19,  longitudinal  and  lateral  distributions  for  10000-keV  down  to 
100-keV  electrons  in  an  infinite  copper  medium  are  presented.  Each  figure  contains  four  results 
corresponding  to  gradually  increasing  pathlengths  of  100,  300,  1000,  and  3000  mfps.  Again,  the 
most  challenging  problem,  or  longitudinal  and  lateral  distributions  for  10000-keV  electrons  in 
copper,  is  presented  first  in  Figure  14  and  Figure  15.  In  all  cases,  the  hybrid  DCS  utilized  is  in 
excellent  agreement  with  the  analog  benchmark.  For  shorter  pathlengths,  both  discrete  models 
oscillate  about  the  analog  benchmark.  The  oscillations  are  an  effect  resulting  from  the 
discreteness  of  the  DCS  model  where  electrons  tend  to  travel  in  preferential  directions.  Even  for 
pathlengths  of  1000  analog  mfps,  the  16-angle  discrete  model  still  oscillates  subtly  about  the 
analog  benchmark,  but  does  not  for  pathlengths  of  3000  analog  mfps. 


Figure  14:  Comparison  of  longitudinal  distributions  for  10000-keV  electrons  after  traveling 
a  distance  of  100  (top  left),  300  (top  right),  1000  (bottom  left),  and  3000  (bottom  right) 
analog  elastic  mfps. 
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Figure  15:  Comparison  of  lateral  distributions  for  10000-keV  electrons  after  traveling  a 
distance  of  100  (top  left),  300  (top  right),  1000  (bottom  left),  and  3000  (bottom  right)  analog 
elastic  mfps. 

As  the  energy  of  the  particle  decreases,  the  effectiveness  of  relaxed  models  (that  is,  models 
preserving  fewer  moments)  is  improved.  For  example,  longitudinal  and  lateral  distributions  for 
1000-keV  electrons  in  copper  are  presented  first  in  Figure  16  and  Figure  17.  In  all  cases,  the 
hybrid  DCS  utilized  is  in  excellent  agreement  with  the  analog  benchmark  and  was  relaxed  from  a 
cutoff  of  0.99  down  to  a  cutoff  of  0.9  for  pathlengths  of  100  and  300  analog  mfps  and  down  to  a 
cutoff  of  0.5  for  pathengths  of  1000  and  3000  analog  mfps.  Moreover,  discrete  artifacts  are  not 
nearly  as  significant  for  1000-keV  electrons  with  exception  of  longitudinal  and  lateral 
distributions  for  pathlengths  of  100  and  300  analog  mfps.  In  fact,  a  4-angle  discrete  model  is 
sufficient  for  resolving  longitudinal  and  lateral  distributions  for  pathlengths  of  1000  and  3000 
analog  mfps. 
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Figure  16:  Comparison  of  longitudinal  distributions  for  1000-keV  electrons  after  traveling 
a  distance  of  100  (top  left),  300  (top  right),  1000  (bottom  left),  and  3000  (bottom  right) 
analog  elastic  mfps. 
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Figure  17:  Comparison  of  lateral  distributions  for  1000-keV  electrons  after  traveling  a 
distance  of  100  (top  left),  300  (top  right),  1000  (bottom  left),  and  3000  (bottom  right)  analog 
elastic  mfps. 
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Finally  for  100-keV  electrons,  the  most  relaxed  models  tested  are  sufficient  for  pathlengths  down 
to  100  analog  mfps  as  seen  in  Figure  18  and  Figure  19.  Although  single-angle  results  were  not 
presented,  a  single-angle  model  is  reasonably  accurate  for  100-keV  electrons  for  pathlengths  of 
100  analog  mfps  and  greater. 


Figure  18:  Comparison  of  longitudinal  distributions  for  100-keV  electrons  after  traveling  a 
distance  of  100  (top  left),  300  (top  right),  1000  (bottom  left),  and  3000  (bottom  right)  analog 
elastic  mfps. 
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Figure  19:  Comparison  of  lateral  distributions  for  100-keV  electrons  after  traveling  a 
distance  of  100  (top  left),  300  (top  right),  1000  (bottom  left),  and  3000  (bottom  right)  analog 
elastic  mfps. 


In  this  section,  longitudinal  and  lateral  results  were  presented  to  demonstrate  the  effectiveness  of 
the  MP  method  when  calculating  quantities  that  are  critical  to  most  CH  methods  and  electron 
transport  methods  in  general.  It  was  shown  that  the  average  longitudinal  and  lateral  displacement 
(typically  used  in  CH  pathlength  correction  algorithms)  is  in  exact  agreement  with  analog  results 
for  a  single-angle  (two  moment-preserving)  model.  In  other  words,  no  additional  pathlength 
correction  algorithm  is  required  in  the  MP  method  because  the  ROP  DCS  are  constructed  such 
that  the  Lewis  moments  are  inherently  preserved.  In  addition,  longitudinal  and  lateral 
distributions  for  100-keV  (less-peaked  scattering)  to  10000-keV  (highly-peaked  scattering)  were 
generated  using  various  ROP  models.  Depending  on  the  problem  at  hand,  ROP  models  requiring 
preservation  of  only  a  few  moments  were  required  for  agreement  with  the  analog  benchmark. 
When  necessary,  additional  moments  were  preserved  to  achieve  analog  level  accuracy,  but  there 
was  never  a  problem  too  extreme  that  the  MP  method  failed  to  resolve  the  longitudinal  or  lateral 
distributions. 

5.3.  1-D  and  2-D  Dose  Calculations 

In  this  section,  1-D  and  2-D  dose  results  are  presented.  The  dose  results  were  generated  using  the 
partial-wave  elastic  scattering  DCS  and  the  Moller  inelastic  scattering  model  for  the  analog 
benchmark  and  the  ROP  DCSs  were,  in  turn,  constructed  from  these  DCSs.  Secondary 
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production  was  not  considered  for  this  section.  The  following  results  show  that  the  MP  method 
can  be  used  to  calculate  dose  accurately  in  both  relatively  isotropic  and  highly  peaked  regimes 
regardless  of  the  form  of  the  analog  model  used  to  construct  the  ROP  DCS.  That  is,  accurate 
models  can  be  constructed  from  both  analytical  DCS  and  tabulated  DCS  data.  In  this  report  the 
emphasis  is  on  tabulated  elastic  DCS  data,  as  applications  of  the  MP  method  to  analytical  DCSs 
has  been  demonstrated  in  the  past  [58,  59].  Efficiency  gains  improve  significantly  with 
increasing  source  energies  without  sacrificing  accuracy.  The  trade-of  between  accuracy  and 
efficiency,  where  it  exist,  ultimately  depends  on  the  application  and  the  level  of  accuracy 
required  by  the  user. 

5.3.1.  One-Dimensional  Depth-Dose  Profiles 

Transversely  integrated  depth-dose  profiles  and  relative  differences  are  presented  for  250-keV 
electrons  in  gold  and  20000-keV  electrons  in  water.  In  Figure  20,  the  results  are  nearly 
indistinguishable  from  the  benchmark.  However,  the  relative  error  plot  in  Figure  20  shows 
disagreements  that  would  otherwise  be  indistinguishable.  In  this  simulation,  disagreement  is 
attributed  to  both  the  elastic  and  inelastic  scattering  models.  Additional  angles  or  use  of  the 
hybrid  DCS  smooths  out  the  overestimation  in  the  first  cell.  Adding  another  energy  point  to  the 
discrete  inelastic  model  smooths  out  the  oscillation  that  begins  near  the  peak  dose.  Though  some 
refinement  was  necessary,  only  small  adjustments  were  required  to  reduce  the  relative 
differences  to  within  1%. 


Figure  20:  Dose  deposition  profile  (left)  and  relative  error  (right)  for  250-keV  electrons  on 
a  gold  slab. 


Efficiency  gains  are  presented  in  Tables  9.5  and  9.6.  For  250-keV  electrons  in  gold  efficiencies 
range  from  about  3  to  50  times  faster  than  analog  depending  on  the  accuracy  of  the  ROP  DCS 
used.  For  20000-keV  electrons,  in  gold  efficiencies  range  from  about  70  to  1800  times  faster 
than  analog,  while  achieving  accuracies  nearly  the  same  as  those  presented  in  Fig.  9. 18b. 
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Table  5:  Efficiency  gains  for  various  discrete  DCSs  when  calculating  dose  due  to  250-kev, 
1000-keV,  and  20000-keV  electrons  on  gold  slabs. 


Particle  Energy 

Reduced  Order  Physics  Model 

1 -Angle 

1 -Energy 

2-Angles 

1 -Energy 

2-Angles 

2-Energies 

4-Angles 

1 -Energy 

4-Angles 

4-Energies 

250-keV 

51 

23 

23 

11 

11 

1000-keV 

164 

72 

71 

31 

31 

20000-kev 

1794 

967 

883 

416 

379 

Table  6:  Efficiency  gains  for  various  hybrid  DCSs  when  calculating  dose  due  to  250-kev, 
1000-keV,  and  20000-keV  electrons  on  gold  slabs. 


Reduced  Order  Physics  Model 


Particle  Energy 

Mo  =  0.5 

1 -Angle 
2-Energies 

Mo  =  0.9 

1 -Angle 

2 -Energies 

Mo  =  0.99 

1 -Angle 
2-Energies 

250-keV 

18 

7 

3 

1000-keV 

60 

22 

6 

20000-kev 

844 

302 

68 

In  Figure  21,  the  results  are  distinguishable  from  the  benchmark  at  20000-keV  because  the 
relaxed  approximations  do  not  capture  large  angle  scatter  or  large  energy  losses.  Once  again, 
small  refinements  to  the  elastic  and  inelastic  ROP  scattering  models  improve  the  accuracy  of  the 
results.  As  seen  in  Figure  21,  model  refinement  through  preservation  of  additional  moments 
reduces  the  relative  differences  to  <1%.  Efficiency  gains  for  the  20000-keV  water  simulation 
range  from  about  1 10  to  1600.  However,  accuracies  within  1%  were  achieved  with  models  that 
were  650  to  700  times  more  efficient  than  analog.  Additional  efficiency  gains  are  presented  in 
Tables  9.7  and  9.8. 


Figure  21:  Dose  deposition  profile  and  relative  error  for  20000-keV  e"  on  H20. 
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Table  7:  Efficiency  gains  for  various  discrete  DCSs  when  calculating  dose  due  to  250-kev, 
1000-keV,  and  20000-keV  electrons  on  water  slabs. 


Particle  Energy 

Reduced  Order  Physics  Model 

1 -Angle 

1 -Energy 

2-Angles 

1 -Energy 

2-Angles 

2-Energies 

4-Angles 

1 -Energy 

4-Angles 

4-Energies 

250-keV 

54 

30 

24 

15 

11 

1000-keV 

147 

83 

61 

39 

27 

20000-kev 

1607 

1122 

709 

612 

293 

Table  8:  Efficiency  gains  for  various  hybrid  DCSs  when  calculating  dose  due  to  250-kev, 
1000-keV,  and  20000-keV  electrons  on  water  slabs. 


Reduced  Order  Physics  Model 


Particle  Energy 

Mo  =  0.5 

1 -Angle 
2-Energies 

Mo  =  0.9 

1 -Angle 

2 -Energies 

Mo  =  0.99 

1 -Angle 
2-Energies 

250-keV 

20 

9 

3 

1000-keV 

52 

25 

7 

20000-kev 

647 

380 

109 

In  addition  to  the  single-material  depth-dose  profiles,  an  interface  problem  is  presented.  In  this 
problem  a  150-keV  pencil  beam  of  electrons  is  normally  incident  on  a  gold-aluminum  slab.  The 
first  0.0004  cm  of  the  slab  is  gold  and  the  remainder  of  the  slab  is  aluminum.  In  Figure  22,  the 
depth-dose  profiles  for  the  analog  benchmark  and  a  single-angle,  single-energy  discrete  model 
are  presented  along  with  the  relative  error  in  several  discrete  models.  The  interface  occurs 
between  the  4th  and  5th  cells;  however,  there  is  no  distinguishable  error  in  Figure  22  resulting 
from  the  interface.  As  previously  noted,  the  MP  method  preserves  transport  mechanics  allowing 
for  exponentially  distributed  collision  sites.  Clearly,  boundary  crossings  are  a  non-issue  for  this 
method  and  no  additional  algorithms  are  required  to  handle  boundary  crossings. 


Figure  22:  Dose  deposition  profile  (left)  and  relative  error  (right)  for  150-keV  electrons  on 
a  gold/aluminum  slab. 
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5.3.2.  Two-Dimensional  Dose  Deposition 

As  indicated  by  the  transversely  integrated  dose  results  in  the  previous  section,  this  method 
provides  excellent  accuracy  and  efficiency  when  radial  spreading  is  not  considered.  However,  in 
this  section  we  present  two-dimensional  dose  deposition  results  that  include  the  impact  of 
radially  spreading.  At  lower  energies  this  method  is  more  effective  at  capturing  radial  spreading. 
We  begin  the  two-dimensional  dose  deposition  results  by  presenting  the  low  energy  simulation. 
The  geometry  setup  for  this  simulation  is  presented  in  Figure  23.  For  the  150-keV  simulation 
109  histories  were  completed.  In  Figure  23,  the  analog  benchmark  for  the  150-keV  simulation  is 
presented.  A  significant  portion  of  the  dose  is  deposited  along  the  beamline  close  to  the  source. 
In  the  gold  region  the  dose  is  deposited  more  rapidly  than  in  the  silicon  region  where  it  is 
apparent  that  the  dose  diffuses  slower. 


The  pencil  beam  is  normally  incident 
at  the  midpoint  of  the  bottom  of  the  block. 
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Figure  23:  Problem  setup  (left)  and  analog  benchmark  (right)  for  150-keV  electrons  on  a 
silicon  cube  with  gold  region. 


In  Figure  24,  relative  differences  for  discrete  models  are  presented  to  demonstrate  the  impact  of 
adjusting  the  number  of  discrete  angles  and  energies,  and  in  turn,  the  number  of  moments 
preserved.  In  Figure  24  (top),  low-order  moment-  preserving  models  are  presented.  Discrete 
artifacts  are  very  distinct  in  Figure  24  (top  left)  for  the  single-angle,  single-energy  model.  By 
increasing  the  number  of  discrete  angles  in  Figure  24  (bottom),  the  discrete  artifacts  are 
mitigated  without  requiring  the  hybrid  model.  However,  there  are  still  some  significant 
differences  in  the  dose  in  some  regions  resulting  from  the  single-energy  model  seen  in  Figure  24 
(top  right  and  bottom  left).  By  including  additional  energies,  the  relative  error  in  all  regions  is 
significantly  reduced  as  seen  in  Figure  24  (bottom  right).  Notice  that  along  the  gold-silicon 
interface  near  the  beam,  the  agreement  is  within  1%,  with  exception  of  the  single-angle  model 
where  the  solution  is  overwhelmed  with  discrete  artifacts. 


Another  artifact  to  point  out  is  the  over/under  estimation  of  the  dose  in  the  first  two  cells  next  to 
the  source.  This  error  persists  as  the  models  are  refined,  but  it  is  actually  an  artifact  of  the  source 
type  and  source  location.  That  is,  the  source  is  a  pencil  beam  that  is  singular  in  space  and  it  is 
directed  at  the  tally  cell  boundary  between.  In  the  next  2-D  result,  the  source  is  no  longer 
singular  in  space  and  the  artifact  is  no  longer  present. 
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PosJUon  (cm)  Position  (on) 

Figure  24:  Relative  error  in  discrete  1-Angle,  1-Energy  model  (top  left),  discrete  2-Angles, 
1-Energy  model  (top  left),  discrete  4- Angles,  1-Energy  model  (top  left),  discrete  4-Angles,  4- 
Energies  model  (top  left)  for  150-keV  electrons  on  silicon/gold  cube. 

Table  9:  Efficiency  gains  for  various  discrete  DCSs  when  calculating  dose  due  to  150-kev 
electrons  on  silicon/gold  cube. 


Reduced  Order  Physics  Model 

Particle  Energy 

1 -Angle 

1 -Energy 

2-Angles 

1 -Energy 

2-Angles 

2 -Energies 

4-Angles 

1 -Energy 

4-Angles 

4-Energies 

250-keV 

51 

29 

29 

16 

16 

In  Figure  25  (left),  the  two-dimensional  problem  setup  is  given  along  with  the  analog 
benchmark  in  Figure  25  (right).  In  this  simulation,  a  beam  of  10000-keV  electrons  are 
transported  with  radius  of  0.02  cm  is  normally  incident  on  a  water  cube  with  a  small  bone  region. 
The  analog  benchmark  is  in  logscale  and  provides  a  sense  of  where  most  of  the  dose  is  deposited. 
That  is,  a  significant  portion  of  the  dose  is  deposited  along  the  beamline  within  fractions  of  a  cm 
to  the  left  and  right  of  the  origin.  The  electrons  with  these  energies  penetrate  deeply  into  the 
medium  as  seen  in  Figure  25  (right). 
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Figure  25:  Problem  setup  (left)  and  analog  benchmark  (right)  for  10000-keV  electrons  on  a 
water  cube  with  bone  region. 

The  following  figure  presents  relative  error  results  corresponding  two  discrete  ROP  DCS 
models.  In  Figure  26  (left),  the  relative  error  in  the  four-angle,  four-energy  discrete  DCS  models 
with  respect  to  the  analog  benchmark  is  presented.  Discrete  artifacts  can  be  seen  clearly  in 
Figure  26  (left).  However,  by  refining  the  model  through  preservation  of  additional  moments 
with  the  addition  of  four  more  discrete  angles,  discrete  artifacts  are  mitigated  and  backscatter  is 
captured  more  accurately  as  seen  in  the  relative  error  result  in  Figure  26  (right).  In  both  the  four- 
angle  and  eight-angle  results,  no  interface  effects  are  present.  It  should  be  noted  that  in  the 
backscatter  is  not  significant  at  10000-keV  and  some  of  the  error  in  the  lower  left  and  right 
comers  is  statistical  in  nature. 


Position  (cm)  Position  (cm) 


Figure  26:  Relative  error  in  discrete  4-Angles,  4-Energies  model  (left)  and  discrete  8- 
Angles,  4-Energies  model  (right)  for  10000-keV  electrons  on  water/bone  cube. 

In  Figure  27,  the  relative  error  between  the  analog  benchmark  and  two  hybrid  models  are 
presented.  In  Figure  27  (left),  a  hybrid  model  with  p*  =  0.9  shows  very  subtle  discrete  artifacts, 
but  otherwise  is  in  good  agreement.  In  Figure  27  (right),  a  hybrid  model  with  p*  =  0.99  does  not 
suffer  from  any  discrete  artifacts  and  the  only  disagreement  is  in  the  lower  left  and  right  comers 
where  again  the  error  is  statistical  in  nature. 
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In  adding  more  discrete  angles,  or  through  use  of  the  hybrid  model,  we  demonstrated  an 
approach  to  mitigate  discrete  artifacts  through  controlling  the  accuracy  of  the  ROP  DCS  models. 
The  only  drawback  to  applying  a  more  accurate  ROP  DCS  model  is  the  loss  of  efficiency  as 
presented  in  Table  9.10.  Here,  the  4-angle  model  is  the  most  efficient  as  expected.  Inclusion  of 
additional  discrete  angles  or  use  of  the  hybrid  model  reduces  the  efficiency  gain  from  roughly 
120  to  between  47  and  94. 


Position  (cm)  Position  (cm) 

Figure  27:  Relative  error  in  hybrid  1-Angle  with  cutoff  of  0.9,  discrete  4-Energies  model 
(left)  and  hybrid  1-Angle  with  cutoff  of  0.99,  discrete  4-Energies  model  (right)  for  10000- 
keV  electrons  on  water/bone  cube. 

Table  10:  Efficiency  gains  for  various  discrete  DCSs  when  calculating  dose  due  to  10000- 
keV  electrons  on  water/bone  cube. 


Reduced  Order  Physics  Model 

Particle  Energy 

4-Angles  8-Angles 

4-Energies  4-Energies 

1 -Angle  1 -Angle 

Mo  =  0.9  fi*0  =  0.99 
4-Energies  4-Energies 

As 

previously 
noted,  it  is 

10000-keV 

121  86 

94  47 

possible  to 

optimize 


such  that  significant  reduction  in  efficiency  is  not  incurred  by  applying  higher-order  models  in 
regions  nearby  the  source  where  the  solution  remains  highly  peaked  and  splaying  lower-order 
models  in  regions  where  the  solution  is  less-peaked.  The  following  results  present  region 
dependent  elastic  ROP  DCSs. 

In  Figure  28,  the  schematic  for  a  region  dependent  discrete  elastic  model  is  presented  along  with 
the  associated  relative  error  from  such  an  approach.  As  seen  in  Figure  28  (left),  default  4-angle, 
4-energy  discrete  model  is  applied  to  all  regions.  The  default  is  then  deactivated  in  the  region 
bounded  by  the  red  dashed  line  and  an  8-angle,  4-energy  discrete  model  is  applied  in  this  region. 
The  associated  relative  error  is  given  in  Figure  28  (right).  The  relative  error  in  Figure  28  (right) 
is  nearly  indistinguishable  from  the  relative  error  in  Figure  28  (right)  where  an  8-angle  discrete 
model  is  applied  everywhere.  The  resulting  gain  in  efficiency  is  106  times  faster  than  the  analog 
simulation  as  opposed  to  86  when  applying  an  8-angle  model  to  all  regions. 
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Figure  28:  Schematic  for  region  dependent  ROP  DCS  (left)  and  the  relative  error  in  dose 
(right)  for  10000-keV  electrons  on  a  water  cube  with  bone  region. 


In  Figure  29,  the  schematic  for  another  region  dependent  discrete  elastic  model  is  presented 
along  with  the  associated  relative  error  from  such  an  approach.  Again,  in  this  problem  a  4-energy 
discrete  inelastic  model  is  used  in  all  regions.  As  seen  in  Figure  29  (left),  an  8-angle  model  is 
applied  in  the  region  where  the  peak  dose  occurs  and  a  single-angle  model  is  applied  in  all  other 
regions.  The  associated  relative  error  is  given  in  Figure  29  (right)  and  is  a  modest  improvement 
over  the  relative  error  in  Figure  29  (right)  as  the  backscatter  is  captured  more  accurately.  Again, 
the  relative  error  in  Figure  29  (right)  is  nearly  indistinguishable  from  the  relative  error  in  Figure 
29  (right)  where  an  8-angle  discrete  model  is  applied  everywhere.  The  resulting  gain  in 
efficiency  is  97  times  faster  than  the  analog  simulation  and  is  reduced  slightly  from  106,  which 
was  efficiency  gain  associated  with  the  previous  region  dependent  models,  but  accuracy  was 
improved. 
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Figure  29:  Schematic  for  region  dependent  ROP  DCS  (left)  and  the  relative  error  in  dose 
(right)  for  10000-keV  electrons  on  a  water  cube  with  bone  region. 

Ultimately,  a  region  dependent  application  of  the  ROP  DCS  models  is  simply  an  exercise  in 
demonstrating  that  accuracy  and  efficiency  can  be  optimized.  In  practice,  a  more  suitable 
approach  would  be  to  develop  an  algorithm  that  determines  the  optimal  DCS  model  as  the 
electron  is  transported.  At  this  point,  it  is  unclear  what  metric  should  is  ideal  for  determining  the 
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optimal  ROP  DCS  model  because  the  parameter  space  is  large.  The  following  problem  is 
illustrates  this  point. 

For  problems  with  singular  boundary  conditions  like  the  previous  problem,  it  is  clear  that  high- 
order  models  are  necessary  nearby  the  source  and  low-order  models  can  be  used  away  from  the 
source  as  the  solution  becomes  less  peaked.  However,  the  same  is  not  necessarily  true  for 
distributed  sources.  For  example,  in  the  following  problem  an  isotropic  point  source  of  2500-keV 
electrons  in  a  gold  cube  1  cm  on  each  face  is  simulated.  The  problem  setup  and  analog  benchmark 
is  given  in  Figure  30  (left).  As  seen  in  Figure  30  (right),  the  dose  is  deposited  uniformly  about  the 
point  source  located  at  the  origin.  It  is  also  of  interest  to  point  out  that  most  of  the  dose  is  deposited 
nearby  the  source. 


0 
0 

0. 

0 
0 

0 

0 
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Figure  30:  Problem  setup  (left)  and  analog  benchmark  (right)  for  an  isotropic  source  of 
2500-keV  electrons  in  a  gold  cube. 

As  seen  in  Figure  31,  relatively  low-order  models  can  used  to  estimate  the  dose  due  to  distributed 
sources.  The  most  efficient  model  tested,  single-angle,  single-energy,  is  about  400  times  more 
efficient  than  the  analog  simulation  and  the  associated  relative  error  is  presented  in  Figure  31 
(left).  By  adding  another  discrete  energy  the  relative  error  improves  as  seen  in  Figure  31  (middle), 
but  the  efficiency  decrease  to  about  300  times  faster  than  the  analog  simulation.  Finally,  the  most 
accurate  model  tested,  2-angles,  2-energies,  is  presented  in  Figure  31  (right).  This  model  provides 
good  agreement,  while  remaining  roughly  180  times  more  efficient  than  the  analog  simulation.  Here 
we  showed  that  accuracy  and  efficiency  is  impacted  by  the  source  configuration.  Again,  region 
dependent  models  could  be  applied  in  this  setting  for  optimization,  but  more  importantly  an  adaptive 
cross-section  algorithm  that  incorporates  source  information  and  solution  information  for  given 
problem  would  improve  the  MP  method. 


1  cm 
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Figure  31:  Relative  error  in  discrete  1-Angle,  1-Energy  model  (left),  discrete  2-Angles,  1- 
Energy  model  (middle)  and  discrete  2-Angles,  2-Energies  model  (right)  for  2500-keV 
electrons  in  gold  cube. 

5.4.  Comparison  with  Experiment 

In  this  section,  an  initial  validation  of  the  Moment-Preservation  method  is  presented.  Given  the 
nature  of  the  method,  the  key  is  to  first  obtain  analog  elastic  and  inelastic  DCSs  that  are  in  good 
agreement  with  the  experimental  benchmarks  of  interest.  Therefore,  the  following  results  are 
really  a  validation  of  the  analog  DCS  models  used  herein.  Initial  validation  results  indicate  that 
the  renormalized  Moller  DCS  (see  Chapter  3)  does  not  give  good  agreement  with  the  Lockwood 
energy  deposition  data  [72].  Use  of  the  Geant4  G4eIonization  class  in  place  of  the  renormalized 
Moller  DCS  improved  agreement.  However,  when  comparing  with  the  Tabata  charge  deposition 
data  [73],  differences  between  the  renormalized  Moller  DCS  and  the  G4eIonization  class  were 
negligible.  Until  further  development  of  the  analog  in-  elastic  model  is  completed,  use  is  made  of 
the  G4eIonization  class  (see  Chapter  8)  for  the  validation  test  under  consideration  in  this  section 
to  be  consistent.  Once  again,  the  analog  elastic  DCS  models  used  are  the  partial-wave  elastic 
DCSs  generated  using  the  ELSEPA  code.  The  first  validation  test  includes  comparisons  to  depth- 
dose  profiles  referred  to  as  the  Lockwood  data  [72].  Next,  comparisons  to  charge  deposition 
experiments  due  to  Tabata  [73]  are  presented.  In  addition  to  the  Lockwood  and  Tabata  data, 
numerous  experimental  benchmarks  are  available  for  validation  test  [74,  75,  76,  77,  78,  79,  80, 
81].  However,  further  validation  of  the  MP  method  remains  as  future  work. 

5.4.1.  Energy  Deposition  Profiles 

One  of  the  most  common  electron  transport  results  in  basic  research  is  energy  deposition,  where 
an  accurate  description  of  particle  transport  is  required  for  different  energies  in  various 
materials.  Below,  energy  deposition  profiles  and  total  energy  de-  position  calculations  are 
compared  with  experimental  results  from  Lockwood  et  al.  [71].  The  Lockwood  data  was 
produced  by  Sandia  National  Laboratories  using  a  sophisticated  calorimetric  technique  for 
measuring  absolute,  high-resolution  electron  energy  deposition  profiles  in  a  wide  range  of  materials. 
The  uncertainty  of  the  data  is  estimated  to  be  from  1.0%  to  2.0%. 
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The  comparisons  cover  low-Z  and  high-Z  materials  including  carbon,  aluminum,  molybdenum,  and 
tantalum  for  pencil  beam  sources  with  energies  of  500-keV  and  1000-keV  and  angles  of  incidence 
of  0°  and  60°.  In  each  simulation,  105  source  particles  were  transported  and  the  energy  deposition 
profiles  are  normalized  to  the  mean  CSDA  range  and  the  depth  variable  is  in  terms  of  a  fraction  of 
the  mean  CSDA  range.  The  primary  objective  of  this  comparison  is  to  validate  the  ROP  models  that 
are  the  subject  of  this  report,  but  also  to  demonstrate  how  well  these  models  perform  with  respect  to 
the  current  state-of-art  physics  models  available  to  Geant4  users.  Both  accuracy  and  efficiency 
results  are  presented,  contrasting  the  MP  method  and  the  default  Geant4  electromagnetic  physics. 

Simulations  were  completed  for  three  different  models  where  the  treatment  of  elastic  scattering 
varies  between  each  model.  These  models  include  an  analog  elastic  scattering  model  given  by  the 
partial-wave  DCS,  a  discrete  single-angle  DCS,  and  the  geant4  elastic  multiple  scattering  model 
referred  to  as  the  Urban  model  or  the  G4UrbanMscModel96  class.  Each  of  the  models  tested  used 
the  same  inelastic  scattering  model  and  bremsstrahlung  model  along  with  the  same  physics  for 
transporting  photons  and  positrons.  The  settings  associated  with  the  aforementioned  physics  are  in 
accordance  with  the  Geant4  standard  electromagnetic  physics  list  option  3,  which  was  found  to  give 
the  best  agreement  with  the  Lockwood  data  [82]  and  enforces  the  strictest  multiple  scattering  step 
limitation  [83]  (this  is  only  relevant  to  the  Urban  model).  In  addition,  the  maximum  step-size  was 
set  to  0.01  mm  for  carbon  and  aluminum  and  0.001  mm  for  molybdenum  and  tantalum  when  using 
the  Urban  model  (step  limits  are  not  required  for  the  analog  model  or  the  discrete  model). 

In  general,  the  energy  deposition  profiles  calculated  using  the  analog  model,  the  ROP  model,  and 
the  geant4  CH  model  exhibit  behavior  similar  to  the  experimental  results.  The  Geant4  physics 
tends  to  estimate  values  higher  than  the  discrete  model  in  the  peak  energy  deposition  region, 
while  the  discrete  model  tends  to  estimate  higher  values  in  the  tails  of  the  energy  deposition 
profile  (see  Figure  32  through  Figure  38).  This  is  true  for  both  nonnal  incidence  and  60°  of- 
nonnal  incidence.  The  total  energy  deposition  tends  to  be  nearly  the  same  for  both  the  MP 
method  and  the  Geant4  physics  models.  For  normally  incident  electrons,  the  agreement  between 
all  models  and  the  experimentally  detennined  total  energy  deposition  is  roughly  1-3%  relative 
difference  for  all  materials  and  energies  with  exception  of  1000-keV  electrons  on  tantalum, 
which  is  between  roughly  4-5%  (see  Table  11).  For  of-nonnal  incidence,  the  agreement  between 
all  models  and  the  experimentally  detennined  total  energy  deposition  is  roughly  2-4%  relative 
difference  (see  Table  12)  for  all  materials  and  energies  with  exception  of  1000-keV  electrons  on 
molybdenum,  which  is  between  roughly  6-7%.  While  this  level  of  agreement  is  generally 
acceptable,  it  is  of  interest  to  develop  an  analog  inelastic  model  that  improves  overall  agreement 
to  within  a  few  percent  relative  error. 

The  timing  results  are  presented  in  Table  13  and  Table  14.  For  normal  incidence,  efficiency 
gains  are  roughly  the  same  for  the  MP  method  and  the  Geant4  physics.  However,  for  of-nonnal 
incidence,  the  MP  method  is  in  all  cases  a  factor  of  two  times  faster  than  the  Geant4  physics. 

For  energy  deposition  calculations,  neither  the  MP  method  nor  the  Geant4  physics 
overwhelmingly  outperform  the  other.  However,  the  MP  method  was  more  efficient  for  the  of- 
nonnal  incidence  simulations.  In  addition,  if  an  analog  model  were  developed  that  provides 
better  agreement  with  experiment,  the  same  level  of  agreement  would  be  anticipated  for  the  MP 
method.  Therefore,  it  is  important  to  identify  an  ideal  analog  inelastic  model  as  the  partial-wave 
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DCSs  are  assumed  to  be  most  accurate  representation  of  elastic  scattering  and  do  not  contribute 
to  any  disagreement  found  herein. 


Furthennore,  no  multi-region  problems  were  validated  in  this  section.  It  is  known  that  material 
interfaces  are  problematic  for  condensed  history  methods  and  there  are  reported  discrepancies 
between  the  Geant4  physics  and  the  Lockwood  data  for  material  interfaces  [82],  Validation  of 
the  MP  method  for  energy  deposition  in  slabs  with  material  interfaces  remains  as  future  work, 
but  it  has  been  shown  in  the  past  that  the  MP  method  does  not  suffer  from  boundary  crossing 
limitations  [58,  59,  60].  Therefore,  no  significant  interface  discrepancies  are  anticipated. 

Table  11:  Total  energy  deposition  comparison  for  500-keV  and  1000-keV  electrons 
normally  incident  on  aluminum,  molybdenum,  and  tantalum  semi-infinite  slabs. 


Energy 

Material 

Model 

Total  Energy 

Relative 

(keV) 

Type 

Type 

(keV) 

Error 

Analog 

472.921 

-0.013 

500 

aluminum 

Discrete  1 -Angle 

471.169 

-0.016 

G4EmStandard 

468.550 

-0.022 

Analog 

958.381 

-0.012 

1000 

aluminum 

Discrete  1 -Angle 

954.665 

-0.016 

G4EmStandard 

953.089 

-0.017 

Analog 

383.612 

0.028 

500 

molybdenum 

Discrete  1 -Angle 

381.167 

0.022 

G4EmStandard 

376.032 

0.008 

Analog 

801.527 

0.029 

1000 

molybdenum 

Discrete  1 -Angle 

795.534 

0.021 

G4EmStandard 

791.940 

0.017 

Analog 

323.144 

0.023 

500 

tantalum 

Discrete  1 -Angle 

316.823 

0.003 

G4EmStandard 

321.674 

0.018 

Analog 

681.921 

0.052 

1000 

tantalum 

Discrete  1 -Angle 

675.622 

0.043 

G4EmStandard 

678.301 

0.047 
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Table  12:  Total  energy  deposition  comparison  for  500-keV  and  1000-keV  electrons  with  60 
degrees  off-normal  incidence  on  aluminum,  molybdenum,  and  tantalum  semi-infinite  slabs. 


Energy 

Material 

Model 

Total  Energy 

Relative 

(keV) 

Type 

Type 

(keV) 

Error 

Analog 

380.501 

-0.027 

500 

aluminum 

Discrete  1 -Angle 

379.573 

-0.029 

G4EmStandard 

375.752 

-0.038 

Analog 

282.940 

0.040 

500 

molybdenum 

Discrete  1 -Angle 

279.703 

0.028 

G4EmStandard 

277.925 

0.022 

Analog 

597.557 

0.073 

1000 

molybdenum 

Discrete  1 -Angle 

593.249 

0.065 

G4EmStandard 

591.644 

0.062 

Analog 

230.202 

0.014 

500 

tantalum 

Discrete  1 -Angle 

226.327 

-0.003 

G4EmStandard 

232.331 

0.023 

Analog 

493.342 

0.045 

1000 

tantalum 

Discrete  1 -Angle 

489.207 

0.036 

G4EmStandard 

495.368 

0.050 
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Table  13:  Timing  results  for  energy  deposition  calculations  for  500-keV  and  1000-keV 
electrons  normally  incident  on  carbon,  aluminum,  molybdenum,  and  tantalum  semi¬ 
infinite  slabs. 


Energy 

Material 

Model 

CPU  time 

Efficiency 

(keV) 

Type 

Type 

(mins) 

Gains 

Analog 

75 

1 

1000 

carbon 

Discrete  1 -Angle 

3 

25 

G4EmStandard 

3 

25 

Analog 

41 

1 

500 

aluminum 

Discrete  1 -Angle 

3 

14 

G4EmStandard 

2 

21 

Analog 

74 

1 

1000 

aluminum 

Discrete  1 -Angle 

4 

19 

G4EmStandard 

3 

25 

Analog 

100 

1 

500 

molybdenum 

Discrete  1 -Angle 

2 

50 

G4EmStandard 

2 

50 

Analog 

151 

1 

1000 

molybdenum 

Discrete  1 -Angle 

3 

50 

G4EmStandard 

5 

30 

Analog 

90 

1 

500 

tantalum 

Discrete  1 -Angle 

2 

45 

G4EmStandard 

1 

90 

Analog 

134 

1 

1000 

tantalum 

Discrete  1 -Angle 

3 

45 

G4EmStandard 

3 

45 
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Table  14:  Timing  results  for  energy  deposition  calculations  for  500-keV  and  1000-keV 
electrons  with  60  degrees  off-normal  incidence  on  aluminum,  molybdenum,  and  tantalum 
semi-infinite  slabs. 


Energy 

Material 

Model 

CPU  time 

Efficiency 

(keV) 

Type 

Type 

(mins) 

Gains 

Analog 

60 

1 

500 

aluminum 

Discrete  1 -Angle 

2 

30 

G4EmStandard 

5 

12 

Analog 

71 

1 

500 

molybdenum 

Discrete  1 -Angle 

1 

71 

G4EmStandard 

2 

35 

Analog 

115 

1 

1000 

molybdenum 

Discrete  1 -Angle 

2 

58 

G4EmStandard 

4 

29 

Analog 

61 

1 

500 

tantalum 

Discrete  1 -Angle 

1 

61 

G4EmStandard 

2 

30 

Analog 

103 

1 

1000 

tantalum 

Discrete  1 -Angle 

1 

103 

G4EmStandard 

2 

52 

Figure  32:  Comparison  with  Lockwood  data  for  1000-keV  electrons  normally  on  carbon 
slab. 
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Figure  33:  Comparison  with  Lockwood  data  for  500-keV  (left)  and  1000-keV  (right) 
electrons  normally  on  aluminum  slab. 


Figure  34:  Comparison  with  Lockwood  data  for  500-keV  (left)  and  1000-keV  (right) 
electrons  normally  on  molybdenum  slab. 


Fraction  of  a  Mean  Range 


Figure  35:  Comparison  with  Lockwood  data  for  500-keV  (left)  and  1000-keV  (right) 
electrons  normally  on  tantalum  slab. 
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Fraction  of  a  Mean  Range 


Figure  36:  Comparison  with  Lockwood  data  for  1000-keV  electrons  with  60  degrees  off- 
normal  incidence  on  aluminum  slab. 


Fraction  of  a  Mean  Range 


Figure  37:  Comparison  with  Lockwood  data  for  500-keV  (left)  and  1000-keV  (right)  with 
60  degrees  off-normal  incidence  on  molybdenum  slab. 
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Figure  38:  Comparison  with  Lockwood  data  for  500-keV  (left)  and  1000-keV  (right) 
electrons  with  60  degrees  off-normal  incidence  on  tantalum  slab. 

5.4.2.  Charge  Deposition  Profiles 

Another  important  electron  transport  result  is  charge  deposition.  Charge  deposition  is  important 
to  understanding  charge  buildup  in  nonconductive  materials.  Below,  charge-deposition  profdes 
are  compared  with  experimental  results  from  Tabata  et  al.  [73].  Among  the  published 
experimental  results  of  charge  deposition  distributions,  those  of  Tabata  et  al.  [74]  cover  the 
widest  regions  of  absorber  atomic  number  (from  4  to  79)  and  incident-electron  energy  (from  4.09 
to  23.5  MeV)  [73]. 

Calculations  of  charge  deposition  distributions  were  performed  for  normally  incident  electron 
pencil  beams  with  energies  of  5000-keV,  10000-keV,  and  20000-keV  on  aluminum  and  gold 
targets  of  thickness  2.5r0,  where  rO  is  the  CSDA  range.  A  total  of  2.4E+05  source  particles  were 
simulated.  The  target  regions  were  divided  into  80  sub-regions  for  scoring  (scoring  regions  are 
not  physical  boundaries)  for  5000-keV  source  particles  and  40  sub-regions  for  10000-keV  and 
20000-keV  source  particles.  Particles  are  tracked  down  to  250-keV. 

Figure  39  through  Figure  44,  present  comparisons  of  charge-deposition  and  dose  distributions 
generated  using  the  partial-wave  elastic  DCS  or  a  discrete  DCS  with  the  default  Geant4  inelastic 
physics  for  electrons  and  the  default  Geant4  physics  for  positrons  and  photons.  The  charge- 
deposition  distributions  are  compared  with  experimental  results  from  Tabata  et  al.  [73].  In 
general  the  correct  behavior  of  the  charge-  deposition  distribution  is  captured  for  all  energies  and 
materials  tested.  In  aluminum,  agreement  with  the  experimental  benchmark  is  satisfactory  with 
exception  of  a  slight  shift  in  the  distribution.  In  gold,  agreement  with  the  experimental  bench¬ 
mark  is  again  satisfactory  with  exception  of  the  5000-keV  results  where  the  charge  deposition  is 
overestimated.  Notice  that  the  analog  solution  and  the  discrete  solution  are  in  excellent 
agreement.  Again,  the  conclusion  is  that  an  analog  model  that  gives  good  agreement  should  be 
identified  and  it  is  assumed  that  the  ROP  models  will  provide  similar  levels  of  agreement. 

Further  investigation  is  required  to  understand  the  disagreements  seen  below  (in  particular,  5000- 
keV  electrons  on  gold). 
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The  dose  results  in  Figure  39  through  Figure  44,  do  show  good  agreement  in  gold  but  do  not  in 
aluminum.  This  result  was  not  anticipated  considering  the  level  of  agreement  in  the  previous 
section  and  suggests  that  future  validations  should  be  avoided  until  this  disagreement  can  be 
identified.  Possible  sources  of  the  disagreement  includes  the  increased  presence  of 
bremsstrahlung  radiation  at  higher  energies,  an  inconsistency  between  the  Geant4  physics 
implementation  and  discrete  model,  or  a  coding  error.  Regardless,  future  work  will  include  the 
use  of  an  inelastic  ROP  model  based  on  an  improved  inelastic  DCS,  so  that  reliance  on  the 
Geant4  physics  is  reduced. 

Table  15:  Timing  results  for  charge  deposition  calculations  for  5000-keV,  10000-keV,  and 
20000-keV  electrons  with  normally  incident  on  aluminum  and  gold  semi-infinite  slabs. 


Energy 

Material 

Model 

CPU  time 

Efficiency 

(keV) 

Type 

Type 

(hours) 

Gains 

Analog 

1.8 

1 

5000 

aluminum 

Discrete  4-Angles 

0.009 

178 

Discrete  8-Angles 

0.02 

80 

Analog 

2.7 

1 

5000 

aluminum 

Discrete  4-Angles 

0.01 

270 

Discrete  8-Angles 

0.03 

90 

Analog 

5.8 

1 

5000 

aluminum 

Discrete  4-Angles 

0.02 

290 

Discrete  8-Angles 

0.05 

116 

Analog 

2.3 

1 

5000 

gold 

Discrete  4-Angles 

0.03 

77 

Discrete  8-Angles 

0.06 

38 

Analog 

4.3 

1 

10000 

gold 

Discrete  4-Angles 

0.06 

72 

Discrete  8-Angles 

0.11 

43 

Analog 

13.8 

1 

20000 

gold 

Discrete  4-Angles 

0.2 

69 

Discrete  8-Angles 

0.3 

46 
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Figure  39:  Charge  deposition  comparison  with  Tabata  data  (left)  and  dose  comparison 
(right)  for  5000-keV  electrons  normally  incident  on  aluminum  slab. 


Figure  40:  Charge  deposition  comparison  with  Tabata  data  (left)  and  dose  comparison 
(right)  for  10000-keV  electrons  normally  incident  on  aluminum  slab. 


Figure  41:  Charge  deposition  comparison  with  Tabata  data  (left)  and  dose  comparison 
(right)  for  20000-keV  electrons  normally  incident  on  aluminum  slab. 
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Figure  42:  Charge  deposition  comparison  with  Tabata  data  (left)  and  dose  comparison 
(right)  for  5000-keV  electrons  normally  incident  on  gold  slab. 


Figure  43:  Charge  deposition  comparison  with  Tabata  data  (left)  and  dose  comparison 
(right)  for  10000-keV  electrons  normally  incident  on  gold  slab. 


Figure  44:  Charge  deposition  comparison  with  Tabata  data  (left)  and  dose  comparison 
(right)  for  20000-keV  electrons  normally  incident  on  gold  slab. 
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5.5.  CEASE  Response  Function  Calculation 


Given  the  experimental  validation  and  the  associated  level  of  confidence  in  the  accuracy  of  the 
analog  DCS  models  and  the  ROP  DCSs  constructed  from  such  models,  the  CEASE  response 
function  calculation  is  revisited.  Here,  it  is  shown  that  the  response  function  can  be  generated 
with  an  efficient  ROP  DCS  model  while  remaining  accurate.  For  this  calculation,  three  models 
were  tested  each  with  the  same  positron,  photon,  and  electron  inelastic  physics.  What  varies  for 
each  of  the  models  is  the  elastic  scattering  physics.  An  analog  model  given  by  the  partial-wave 
elastic  DCS  is  tested  and  a  single-angle  discrete  model  based  on  the  partial-wave  elastic  DCS  is 
compared.  In  addition,  the  response  function  is  also  calculated  using  the  default  Geant4 
electromagnetic  physics  option  3.  In  Figure  45,  response  functions  generated  using  each 
different  model  are  compared.  The  analog  and  the  discrete  model  are  in  excellent  agreement, 
while  the  Geant4  physics  (that  is,  class  II  CH  with  out  a  user  applied  step  limitation)  shows 
significant  disagreement  at  higher  energies.  Assuming  the  analog  model  is  the  most  accurate,  the 
Geant4  physics  tends  to  under  predict  the  detector  response  at  higher  energies.  This  could  be  an 
effect  of  the  collimator,  as  it  is  known  that  the  Geant4  physics,  without  step  limitation,  tends  to 
overestimate  energy  deposition  in  high-Z  materials.  Therefore,  it  is  possible  that  electrons  do  not 
fully  penetrate  the  collimator  because  they  lose  too  much  energy  in  the  tungsten  collimator.  See 
Figure  46  for  an  example  of  trajectories  penetrating  the  collimator.  The  collimator  is  in  green; 
the  detectors  are  in  blue;  and  the  electron  tracks  are  in  red. 

The  efficiency  gains  associated  with  the  discrete  DCS  are  outstanding.  To  complete  an  analog 
simulation  of  1 1 16  runs  for  18  source  angles  and  62  source  energies  with  1000  source  particles 
each,  it  requires  roughly  46  processor-days,  while  for  the  discrete  model  it  only  takes  roughly 
five  processor-days.  That  is  an  order  of  magnitude  efficiency  gain.  The  CH  runs  were  roughly  10 
times  faster  than  the  discrete  model  or  0.5  processor-days. 


Figure  45:  Comparison  of  response  functions  generated  using  the  analog  model,  the 
discrete  model,  and  the  Geant4  multiple-scattering  model. 
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Figure  46:  Electrons  traversing  the  CEASE  telescope. 

5.6.  Detector  Modeling  and  Simulation 

The  following  section  present  related  work  including  the  modeling  of  detectors  associated  with 
the  AFRL  Demonstration  and  Science  Experiment  (DsX)  detectors  for  use  in  response  function 
determination.  Modeling  and  analysis  was  completed  for  the  Compact  Environmental  Anomaly 
Sensor  (CEASE),  the  High  Energy  Proton  Spectrometer  (HEPS),  and  the  Low  Energy 
Electrostatic  Analyzer  (LEESA).  Each  detector  was  modeled  using  MCNPX  (Monte  Carlo  N- 
Particle  Extended)  particle  transport  software.  The  detector  models  were  completed  by 
constructing  an  MCNPx  input  file  with  detail  of  the  detector  consistent  with  the  available 
information  (e.g.  detector  drawings,  and  so  on). 

Given  the  highest  possible  fidelity  detector  models,  electron  and  proton  simulations  were 
completed  for  a  wide  range  of  particle  energies  and  angles  of  incidence.  The  output  or  result 
from  the  individual  simulations  was  detailed  particle-by-particle  energy  deposition  for  each 
region  of  interest  (that  is,  sensitive  detector  regions).  In  the  CEASE  detector  the  sensitive  regions 
included  two  silicon  PIPS  detectors;  whereas,  the  HEPS  detector  contains  three  silicon  PIPS 
detectors  and  three  plastic  scintillators,  or  six  sensitive  regions. 

With  detailed  energy  deposition  information  for  the  relevant  source  particle  energies  and  angles, 
and  detector  logic,  the  response  function  for  various  detector  channels  and  or  the  total  response 
function  is  determined.  To  do  so,  a  data  processing  algorithm  is  required  where  the  details  of 
such  an  algorithm  are  given  in  [1].  Ultimately,  the  response  functions  generated  using  models 
built  by  UNM  and  the  UNM  algorithm  for  data  processing  were  satisfactory.  That  is,  good 
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agreement  was  found  between  the  UNM  generated  response  functions  for  CEASE  and  HEPS  and 
relevant  response  functions  from  Brautigan  et  al.  [1], 

The  focus  for  LEESA  was  on  familiarization  with  the  detector  system  and  high  fidelity  modeling 
of  the  detector,  rather  than  conducting  analysis.  Familiarization  with  the  detector  and  an 
extensive  model  was  completed.  Once  the  model  was  completed,  emphasis  was  shifted  to  the 
completing  the  research  under  consideration  for  the  grant. 

6.  CONCLUSIONS 

It  is  of  interest  to  develop  an  alternative  to  the  CH  method  that  is  free  of  the  limitations  inherent  to 
CH.  The  subject  of  this  report,  or  the  Moment-Preserving  method,  is  such  an  alternative  and 
therefore,  the  accuracy  and  efficiency  of  this  method  must  be  demonstrated  and  contrasted  (to  some 
degree)  with  the  CH  method.  To  do  so,  the  theoretical  development  of  the  method  was  discussed  in 
great  detail,  emphasizing  how  elements  of  accuracy  and  efficiency  are  inherent  to  the  method, 
while  also  providing  an  exhaustive  numerical  demonstration  including  validation  of  the  method. 

Through  theoretical  development,  it  was  shown  that  in  the  Moment-Preserving  method  a  reduced 
order  physics  (ROP)  transport  equation  is  formed  by  replacing  the  analog  DCS  with  ROP  DCSs 
that  are  less  peaked  with  longer  mfps.  By  doing  so,  solution  to  the  ROP  transport  equation  using  a 
single-event  Monte  Carlo  method  is  computationally  efficient  with  respect  to  analog  Monte  Carlo. 
To  simply  replace  the  analog  DCSs  with  another  more  ideal  DCS  requires  a  process  for 
constructing  such  a  DCS  that  guarantees  accurate  results.  This  is  achieved  by  applying  a  theory 
ubiquitous  in  electrons  transport  methods,  or  Lewis  theory,  that  relates  moments  of  the  analog  DCS 
to  moments  of  the  solution.  By  recognizing  the  importance  of  Lewis  theory,  a  process  was 
developed  for  constructing  ROP  DCSs  by  applying  a  moment-  preservation  constraint,  where  both 
elastic  and  inelastic  ROP  DCSs  are  constructed  such  that  they  preserve  some  finite  number  of 
moments  of  the  analog  DCS  exactly  and  higher-order  moments  are  approximated  in  terms  of  the 
lower-order  moments.  The  resulting  ROP  DCSs  preserve  key  physical  moments  like  the  mean 
scattering  cosine,  the  transport  cross-section,  the  stopping  power,  and  the  straggling  coefficient, 
along  with  any  other  user  specified  higher-order  moments.  Furthermore,  by  constructing  the  ROP 
DCSs  such  that  one  point  is  required  to  coincide  with  the  nearly-singular  point  of  the  analog  DCS, 
the  method  takes  advantage  of  a  convenient  cancellation  of  in-scatter  and  out-scatter  due  to  these 
nearly-singular  points  about  zero  changes-  in-state  resulting  in  a  reduction  of  the  total  cross- 
section.  Between  the  moment  preserving  constraint  and  the  cancellation  of  the  nearly-singular 
points,  moment-  preserving  ROP  DCSs  that  are  less-peaked  with  a  longer  mfp  (up  to  four  orders  of 
magnitude  longer  than  the  analog  mfp)  can  be  generated  and  accurate  solutions  to  the 
corresponding  ROP  transport  equation  can  be  obtained  efficiently.  These  findings,  though 
satisfactory,  only  partially  satisfy  the  objective  of  the  research,  which  is  to  develop  and  demonstrate 
an  alternative  for  electron  transport. 

To  completely  satisfy  the  objectives,  a  numerical  demonstration  of  the  Moment-Preserving  method 
was  presented.  The  results  associated  with  the  numerical  demonstration,  served  to  extend  and 
modernize  the  significant  body  of  work  completed  by  Franke  and  Prinja  in  2005  [58].  By  extending 
their  study,  several  key  features  of  this  method  were  addressed  including:  systematically 
controllable  accuracy,  efficiency,  mathematical  robustness,  versatility  through  the  independence  of 
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the  method  from  the  form  of  the  analog  DCS,  flexibility  of  the  method  through  usage  of  both  the 
discrete  and  hybrid  DCS,  and  simplicity  easing  implementation  in  transport  codes.  The  accuracy  and 
efficiency  of  the  method  was  demonstrated  through  calculation  of  both  differential  and  integral 
quantities  in  both  highly-peaked  scattering  regimes  and  less-  peaked  scattering  regimes  for  a  wide 
variety  of  target  materials  and  source  energies.  Results  including  reflected  and  transmitted  angular 
distributions  and  energy-loss  spectra  in  thin  slabs,  longitudinal  and  lateral  distributions  in  infinite 
media,  dose  deposition  in  1-D  and  2-D  slabs,  and  charge  deposition  were  presented.  For  the 
theoretical  problems  where  the  method  is  compared  to  an  analog  benchmark,  analog  level 
accuracies  were  achieved  with  efficiency  gains  up  to  three  orders  of  magnitude  greater  than  analog 
level  efficiencies.  For  the  validation  results,  the  accuracies  and  efficiencies  were  similar  to  the 
default  Geant4  electromagnetic  physics  with  a  factor  of  two  improvement  in  efficiency  for  of- 
nonnal  incidence  source  problems.  Although  the  Moment-Preserving  method  was  not  a  dramatic 
improvement  over  the  default  Geant4  electromagnetic  physics,  this  was  a  first  attempt  at  validating 
the  method  and  much  work  remains  in  identifying  an  analog  DCS  model  for  inelastic  scattering. 
Improvements  in  accuracy  beyond  what  was  reported  herein  are  anticipated  with  such  an  inelastic 
DCS  model.  As  it  stands,  the  Moment-Preserving  is  a  suitable  alternative;  however,  with  these 
improvements  one  can  expect  that  the  Moment-Preserving  method  will  provide  a  clear  advantage 
over  the  default  Geant4  electromagnetic  physics. 

In  its  current  state  with  regards  to  accuracy  and  efficiency,  the  Moment-Preserving  method  is  a 
strong  alternative  to  CH  methods.  Under  conditions  where  CH  methods  and  the  Moment-Preserving 
method  provide  identical  levels  of  accuracy  and  efficiency,  the  Moment-Preserving  method  has  the 
added  advantage  of  versatility  and  simplicity.  That  is,  no  changes  to  the  source  code  or  the 
algorithm  are  required  to  make  significant  changes  to  the  underlying  physics.  The  moment¬ 
preserving  algorithms  are  completely  independent  of  the  form  of  the  analog  DCS  that  ultimately 
drives  the  accuracy  of  this  method.  Therefore,  if  an  improved  elastic  or  inelastic  DCS  is 
developed,  one  must  simply  generate  an  ROP  DCS  library  corresponding  to  the  improved  models. 

In  addition,  initial  implementation  of  the  Moment-Preserving  method  is  straightforward;  especially, 
in  transport  codes  with  pre-existing  single-  scatter  algorithms.  As  a  result  of  the  mathematical 
robustness  of  the  method,  no  additional  algorithms  are  required  beyond  what  is  typically  used  for  a 
single-scatter  algorithm  that  uses  tabulated  DCS  data.  The  algorithm  simply  requires  methods  for 
table  look-ups  to  obtain  the  total-cross  section,  methods  for  sampling  DCS  data  without 
interpolation,  and  methods  for  processing  the  ROP  DCS.  In  Geant4,  this  required  utilization  of  the 
existing  data  classes  and  base  classes  for  physics  models  and  processes.  Any  maintenance 
associated  with  Moment-Preserving  method  would  be  a  result  of  changes  to  the  base  classes  from 
which  the  Moment-Preserving  method  physics  inherit. 

Given  the  accuracy,  efficiency,  mathematical  robustness,  versatility,  flexibility,  and  simplicity  of 
the  Moment-Preserving  method,  this  method  provides  a  clear  alternative  to  the  prevailing  electron 
transport  method  -  Condensed  History.  This  work  establishes  a  modern  basis  from  which  further 
testing  of  this  method  will  be  completed. 
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